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Endocrine  therapy  is  often  the  least  toxic  and  most  effective  treatment  for  hormone  receptor  positive  invasive 
breast  cancer.  Such  therapy  includes  antiestrogens  (tamoxifen,  fulvestrant)  and  aromatase  inhibitors 
(anastrozole,  letrozole,  exemestane).  Tamoxifen  (TAM)  increases  disease  free  and  overall  survival  in  the 
adjuvant  setting,  reduces  the  incidence  of  estrogen  receptor  positive  disease  (ER+;  unless  otherwise  noted 
ER=ERa)  in  high-risk  women,  and  reduces  the  rate  of  bone  loss  secondary  to  osteoporosis  in  postmenopausal 
women  [1,2].  Aromatase  inhibitors  are  effective  only  in  the  absence  of  functioning  ovaries  -  TAM  can  be  used 
regardless  of  menopausal  status.  Recent  studies  suggest  that  anastrozole  may  be  superior  to  TAM  in  the 
adjuvant  treatment  of  postmenopausal  women  with  ER+  breast  cancer;  other  studies  report  higher  overall 
response  rates  with  letrozole  (LET)  vs.  TAM  as  first  line  therapy  in  the  metastatic  setting.  Thus,  a  recent 
controversy  in  the  management  of  patients  with  ER+  disease  is  whether  an  aromatase  inhibitor  or  TAM  should 
be  given  as  first  line  endocrine  therapy  [3-9]. 

In  this  Clinical  Translational  Research  award,  we  will  build  classifiers  that  accurately  separate  antiestrogen 
sensitive  from  antiestrogen  resistant  breast  tumors  and  begin  to  assist  in  the  direction  of  specific  endocrine 
treatments  (antiestrogen  vs.  aromatase  inhibitor)  to  individual  patients.  We  hypothesize  that  endocrine 
responsiveness  is  affected  by  a  gene  network,  rather  than  the  activity  of  only  one  or  two  genes  or  signaling 
pathways  [10-12],  Since  the  key  components  of  such  a  network  are  unknown,  we  must  study  10,000s  of  genes. 
We  will  use  Affymetrix  GeneChips.  We  will  not  identify  mutational  events,  the  presence  of  mRNA  splice 
variants,  or  post-translational  protein  modifications.  However,  these  factors  have  major  effects  on  the 
transcriptome  and  their  "footprints"  should  be  identified  by  expression  microarrays. 


Body 

Overview:  We  will  build  classifiers  that  separate  antiestrogen  sensitive  from  antiestrogen  resistant  breast 
tumors  and  begin  to  assist  in  the  direction  of  specific  endocrine  treatments  (antiestrogen  vs.  aromatase  inhibitor) 
to  individual  patients.  To  achieve  this  goal,  and  consistent  with  a  CTR  award,  we  will  complete  a  4-year, 
prospective,  neoadjuvant  study  with  Letrozole  (LET)  or  TAM  as  the  only  systemic  therapy.  We  will  obtain 
molecular  profiles  from  Affymetrix  GeneChips  and  further  develop  and  apply  our  innovative  bioinformatic  and 
biostatistic  methods  to  explore  these  high  dimensional  data  sets  and  build/validate  new  classifiers.  A  more 
accurate  predictor  of  endocrine  responsiveness  would  have  widespread  clinical  use,  allowing  women  and 
physicians  to  make  more  individualized  and  appropriate  treatment  decisions.  For  example,  patients  with  tumors 
predicted  to  be  resistant  to  antiestrogens  and/or  aromatase  inhibitors  would  be  strong  candidates  for  an  early 
intervention  with  cytotoxic  chemotherapy. 

In  most  predictive/prognostic  marker  studies  investigators  focus  on  a  single  factor  and  whether  they  obtain  a  p- 
value  that  reaches  conventional  statistical  significance.  Our  approach  is  different  because  we  will  determine 
whether  we  can  find  joint  gene  subsets  that  can  separate  patients  into  sufficiently  distinct  groups  that  should 
differ  in  their  treatment.  We  will  (1)  analyze  >33,000  genes  on  retrospective  and  prospective  material,  (2)  apply 
new  biostatistical  and  bioinformatic  methods  to  identify  ~40  potentially  informative  "biomarkers,"  (3)  build 
neural  network  and  biostatistical  model  classifiers,  (4)  evaluate  the  joint  discriminant  power  of  selected  genes 
concurrently  rather  than  as  single  biomarkers,  (5)  focus  on  prediction  for  individual  patients  where  the 
assessment  of  a  p- value  is  less  important  than  the  classification  rate  of  our  predictors,  (6)  validate  the  classifiers 
in  independent  data  sets,  and  (7)  explore  the  ability  of  predictors  to  refine  the  targeting  of  specific  endocrine 
therapies. 

Evidence  has  begun  to  accumulate  suggesting  that  an  aromatase  inhibitor  might  be  a  more  effective  first  line 
endocrine  therapy  for  some  breast  cancer  patients  than  the  current  standard  of  care  (Tamoxifen).  These  data 
have  generated  considerable  interest  and  controversy,  in  part  because  unlike  TAM,  there  are  no  long  term 
studies  with  aromatase  inhibitors  where  definitive  survival  data  are  available.  Our  study  could  provide  new  and 
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innovative  insights  into  how  to  approach  the  more  effective  targeting  of  specific  endocrine  therapies  to 
individual  patients. 


Specific  Aims 

We  will  complete  two  clinical  studies  and  collect  gene  expression  profiles  from  which  to  build  predictors  of 
endocrine  responsiveness.  Predictors  will  be  built  in  Specific  Aim  2  and  validated  in  Specific  Aim  3. 

Aim  1 :  Clinical  Studies  -  Clinical  Study-1  (retrospective)  is  of  pretreatment,  single,  frozen  samples  where  we 
will  compare  the  molecular  profiles  of  tumors  that  recurred  on  TAM  with  those  of  tumors  that  did  not  recur. 
Each  resistant  sample  is  matched  with  a  TAM  sensitive  sample  by  age,  stage,  and  duration  of  follow-up.  We 
also  have  further,  single  (unmatched),  frozen  samples  from  patients  already  progressing  on  TAM.  Clinical 
Studv-2  is  a  prospective  study  of  breast  tumor  samples  from  patients  treated  with  neoadjuvant  TAM  or  LET. 

Aim  2:  We  will  develop  and  apply  novel  bioinformatics  and  biostatistics  to  discover  gene  subsets  that  define  the 
molecular  differences  between  endocrine  sensitive  and  resistant  breast  tumors.  These  genes  will  be  used,  in 
combination  with  established  predictive/prognostic  factors,  e.g.,  ER,  PgR,  stage,  to  build  innovative  classifiers 
that  can  better  predict  an  individual  tumor’s  endocrine  responsiveness. 

Aim  3:  We  will  test,  optimize,  and  validate  the  performance  of  the  classifiers  from  Aim  2  in  retrospective 
studies  of  human  breast  tumors.  We  will  measure  each  gene  individually  by  IHC,  in  situ  RNA  hybridization 
(ISH),  or  real  time  PCR  (RT-PCR). 


Key  Research  Accomplishments 

As  noted  in  previous  reports,  progress  on  the  clinical  goals  for  this  award  was  greatly  delayed  because  of  the 
time  taken  to  obtain  DOD  approval  of  our  preexisting  institutionally  approved  IRBs  at  Georgetown  University 
and  at  the  University  of  Edinburgh.  All  institutionally  approved  protocols  and  requested  material  were 
submitted  to  the  DOD  in  July  2004;  additional  information  was  requested  by  the  DOD  several  months  later  and 
submitted  in  November  2004.  We  did  not  receive  final  approval  to  proceed  with  the  clinical  studies  until  March 
2005.  Much  of  this  delay  seems  to  have  been  entirely  unavoidable  (see  prior  reports).  Clearly,  this  has  likely  left 
us  behind  schedule  in  recruitment  to  the  prospective  studies.  We  have  made  significant  strides  in  our 
development  of  new  analytical  procedures,  and  we  have  been  successful  in  using  our  emerging  data  to  support 
additional  applications  for  funding.  Publications  supported  since  the  commencement  of  this  award  are  listed 
under  “Reportable  Outcomes”;  these  constitute  some  of  our  major  accomplishments  in  the  past  year.  These  and 
other  key  research  accomplishments  are  presented  below. 


Progress  on  our  Statement  of  Work  (SOW  from  the  original  application) 

•  TASK  1.  Array  breast  tumor  samples  from  Clinical  Studies  1  (retrospective)  and  2  (prospective) 

We  have  now  received  a  total  of  481  breast  specimens  from  our  collaborators  at  the  University  of  Edinburgh. 
These  specimens  arrived  at  different  times  and  were  initially  banked  so  that  they  could  be  processed  in  the  most 
effective  and  logical  manner.  These  specimens  represent  a  mix  of  the  initial  prospective  and  most  of  the 
retrospective  specimens.  We  have  processed  all  of  these  specimens  as  frozen  sections  and  all  481  have  now 
been  fully  analyzed  and  annotated  by  the  study  pathologist.  We  have  successfully  extracted  total  RNA  from  357 
specimens,  and  labeled  169  for  analysis.  We  have  also  completed  the  hybridization  and  assessment  of 
microarray  data  quality  control  on  102  breast  cancer  specimens. 

We  requested  that  the  specimens  be  sent  independent  of  the  clinical  information,  so  that  we  could 
adequately  and  appropriately  randomize  the  RNA  preparation,  labeling  and  hybridization  and  minimize  any 
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operator-induced  or  technology-induced  bias.  All  specimens  were  processed  using  our  standard  operating 
procedures;  each  manipulation  being  performed  by  the  same  individual  to  further  reduced  inter-operator 
variability.  Details  of  the  methods,  quality  control  measures  and  general  experimental  approaches  have  been 
described  in  detail  in  earlier  annual  reports. 

In  general,  we  have  found  the  material  from  our  collaborators  to  be  of  high  quality.  Over  80%  of  the 
specimens  have  yielded  RNA  with  an  estimated  RIN  value  of  >5.0,  which  we  have  found  produces  good  quality 
expression  data.  Some  of  the  specimens  with  slightly  lower  RIN  values  also  will  be  useful  but  we  will  first 
label,  hybridize  and  array  the  best  quality  specimens.  While  not  all  cases  generate  adequate  RNA,  and  a  small 
number  of  tissues  have  very  little  neoplastic  tissue  (most  likely  the  best  responders),  thee  data  suggest  that  we 
should  be  able  to  address  our  central  hypotheses  as  initially  anticipated. 

We  have  also  found  these  data  to  be  particularly  useful  in  supporting  other  studies  that  are  ongoing  in  the 
laboratory.  For  example,  these  data  have  been  used  to  support  two  pending  R01  applications  on  genes  we 
identified  and  described  in  the  preliminary  data  for  this  application.  We  have  also  used  these  data  to  provide 
preliminary  data  on  gene  expression  values  that  have  led  to  our  colleagues  initiating  other  studies  directed  at 
developing  therapeutic  strategies  to  target  individual  genes  we  have  identified  from  within  this  data  set  or  from 
other  sources. 


•  TASK  2.  Store,  process,  and  train/optimize  classifiers  from  gene  expression  microarray  data  (modified 
to  reflect  our  adoption  of  caArray) 

As  noted  in  our  previous  reports,  we  continue  to  make  significant  progress  on  addressing  this  task,  largely  as  a 
consequence  of  our  involvement  in  the  National  Cancer  Institute  Center  for  Bioinformatics  (NCICB)  led  caBIG 
project.  The  PI  (Dr.  Clarke)  leads  the  Lombardi  Comprehensive  Cancer  Center’s  caBIG  team  and  we  have  been 
actively  involved  in  the  development  of  caArray  (NCICB ’s  grid-enabled,  MIAME  compliant,  microarray 
database).  Indeed,  we  will  host  a  key  caBIG  face-to-face  meeting  at  our  institute  that  will  represent  a  joint 
meeting  between  the  Architecture  and  Vocabulary  and  Common  Data  Elements  Workspaces  members. 

With  respect  to  the  further  development  and  optimization  of  data  analysis  algorithms,  we  have  begun  to  develop 
novel  approaches  for  network  analysis  that  will  allow  us  to  further  identify  potentially  novel  targets  for 
therapeutic  intervention  in  subsequent  studies.  We  anticipated  obtaining  such  information  in  our  original 
application,  and  we  have  found  approaching  this  goal  to  be  realistic  in  a  much  shorter  time  frame  than  initially 
expected. 

We  also  continue  to  improve  our  existing  algorithms  and  have  recently  submitted  for  publication  a  short 
communication  on  the  implementation  and  uses  of  our  VISDA  algorithms  (described  in  the  initial  application). 
As  described  in  last  year’s  report,  we  see  the  potential  to  obtain  novel  mechanistic  insights  as  a  significant 
advantage  to  our  ongoing  studies.  We  will  provide  additional  information  in  this  regard  in  subsequent  reports; 
relevant  publications  in  this  area  are  included  below  in  the  section  “Reportable  Outcomes.” 


•  TASK  3.  Retrain/reoptimize  classifiers  using  IHC  data  from  Series  1  (Archival  Tissues)  and  Series  2 
(Scottish  Adjuvant  TAM  Trial)  for  Validation 

To  perform  this  task  we  will  obtain  clinical  information  and  breast  tumor  samples  from  University  of  Edinburgh 
(formalin  fixed/paraffin  embedded).  We  will  rank  and  prioritize  selected  joint  genes  from  RNA  classifier  built 
and  optimized  in  TASK  2  (above)  and  retrain/reoptimize  the  initial  neural  network  IHC  classifier  (MLP). 
Finally,  we  will  validate  IHC  classifier  on  independent  data  sets  (data  sets  not  used  to  build  and  train  the  MLP 
classifiers). 

We  remain  unable  to  move  this  task  substantially  forward  on  the  timeframe  as  initially  proposed  because  of  the 
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delays  in  getting  approval  to  work  with  the  clinical  specimens.  It  is  addressing  this  aim  in  detail  that  will  be  the 
greatest  consideration  in  assessing  the  need  to  take  the  reviewer’s  advice  from  their  assessment  of  last  year’s 
report  that  a  request  for  a  no-cost  extension  may  be  advisable  (see  also  “Conclusions”  below). 


Reportable  Outcomes 


Papers  and  Meeting  Reports* 

Updates  (cited  as  “in  press”  in  the  last  report  and  now  in  print) 

•  Ressom,  H.W.,  Zhang,  Y.,  Xuan,  J.,  Wang,  Y.  &  Clarke,  R.  “Inference  of  gene  regulatory  networks 
from  time  course  gene  expression  data  using  neural  networks  and  swarm  intelligence.”  IEEE  Symp 
Compl  Intel  Bioinf  Comput  Biol,  435-442,  2006. 


New  Publications  (for  the  present  reporting  period) 

•  Wang,  L.H.,  Yang,  X.Y.,  Zhang,  X.,  An,  P.,  Kim,  H.-J.,  Huang,  J.,  Clarke,  R.,  Osborne,  C.K.,  Inman, 
J.K.,  Appella,  E.  &  Farrar,  W.L.  “Disruption  of  estrogen  receptor  DNA-binding  domain  and  related 
intramolecular  communication  restores  tamoxifen  sensitivity  in  resistant  breast  cancer.”  Cancer  Cell,  10: 
487-499,  2006. 

•  Bouker,  K.B.,  Skaar,  T.C.,  Harburger,  D.S.,  Riggins,  R.B.,  Fernandez,  D.R.,  Zwart,  A.  &  Clarke,  R. 
“The  A4396G  polymorphism  in  interferon  regulatory  factor- 1  is  frequently  expressed  in  breast  cancer.” 
Cancer  Genet  Cytogenet,  175:  61-4,  2007. 

•  Wong,  L.-J.C.,  Dai,  P.,  Lu,  J.-F.,  Lou,  M.A.,  Clarke,  R.  &  Nazarov,  N.  “AIB1  gene  amplification  and 
the  instability  of  polyQ  encoding  sequence  in  breast  cancer  cell  lines.”  BMC  Cancer  6:  111-132,  2006. 

•  Kuske,  B.,  Naughton,  C.,  Moore,  K.,  MacLeod,  K.G.,  Miller,  W.R.,  Clarke,  R.  Langdon,  S.P.  & 
Cameron,  D.A.  “Endocrine  therapy  resistance  can  be  associated  with  high  estrogen  receptor  alpha  (ERa) 
expression  and  reduced  ERa  phosphorylation  in  breast  cancer  models.”  Endocr  Related  Cancer,  13: 
1121-1133,2006. 

•  Gong,  T.,  Xuan,  J.,  Zhu,  J.,  Li,  H.,  Clarke,  R„  Hoffman,  E.  &  Wang,  Y.  “Composite  gene  module 
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2006 
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recurrent  neural  networks  and  swarm  intelligence,"  Proc  28th  IEEE  EMBS  Inti  Conf,  pp.  4241-4244, 
2006. 

*We  include  in  the  appendix  reprints  of  those  papers  that  are  already  published  adn  for  which  we  have  proofs  or 
reprints.  We  do  not  list  here  or  include  in  the  appendices  any  published  abstracts,  but  can  do  so  if  requested. 


Page  7 


Award  Number:  W81XWH-04-1-0570 _ PI:  Robert  Clarke,  Ph.D.,  D.Sc. 

Several  other  manuscripts  related  to  our  bioinformatic  methods  also  are  submitted  and  in  preparation  -  these 
will  be  cited  reported  in  the  next  report.  Please  note  that  the  papers  published  in  the  engineering  literature  are 
different  from  most  conference  proceedings  in  the  biomedical  literature.  These  are  not  abstracts  but  fully  peer- 
reviewed  publications  comparable  to  short  communications  in  biomedical  journals. 

Comment  on  Subcontracts:  Please  also  note  that  the  majority  of  our  publications  here  and  in  prior  years  include 
coauthors  from  one  or  both  of  our  subcontracts.  Thus,  our  program  is  working  very  effectively  and 
collaboratively,  this  should  further  be  apparent  in  the  development  of  new  informatics  methods  (Virginia 
Polytechnic  and  State  University  subcontract)  and  the  large  number  of  high  quality  breast  tumor  specimens  we 
have  obtained  from  the  University  of  Edinburgh. 

Conclusions 

We  have  made  good  progress  on  the  research  infrastructure  goals  and  in  the  development  or  optimization  of  the 
methods  needed  for  data  analysis.  We  also  have  completed  and  published  all  of  the  data  presented  as 
preliminary  data  in  the  initial  application.  The  clinical  studies  were  held  up  by  an  unexpectedly  long  delay  in 
obtaining  final  approval  for  our  existing  protocols.  As  noted  by  the  reviewer  of  last  year’s  annual  report,  this 
delay  will  adversely  affect  the  prospective  study  and  this  reviewer  indicated  that  a  request  for  a  one-year  no  cost 
extension  might  be  required.  We  concur  that  this  may  yet  be  required  and  will  evaluate  the  need  for  this  over 
the  next  six  months.  If  deemed  necessary,  we  will  apply  for  such  an  extension  in  writing  before  the  end  of  the 
original  funding  period,  this  should  ensure  that  we  remain  in  compliance  with  USAMRMC  guidelines,  maintain 
continuity  of  the  project  and  successfully  complete  our  studies. 
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Abstract —  We  present  a  novel  algorithm  that  combines  a 
recurrent  neural  network  (RNN)  and  two  swarm  intelligence  (SI) 
methods  to  infer  a  gene  regulatory  network  (GRN)  from  time 
course  gene  expression  data.  The  algorithm  uses  ant  colony 
optimization  (ACO)  to  identify  the  optimal  architecture  of  an 
RNN,  while  the  weights  of  the  RNN  are  optimized  using  particle 
swarm  optimization  (PSO).  Our  goal  is  to  construct  an  RNN 
whose  response  mimics  gene  expression  data  generated  by  time 
course  DNA  microarray  experiments.  We  observed  promising 
results  in  applying  the  proposed  hybrid  SI-RNN  algorithm  to 
infer  networks  of  interaction  from  simulated  and  real-world  gene 
expression  data. 

I  INTRODUCTION 

Gene  regulatory  network  (GRN)  is  a  model  of  a  network 
that  describes  the  relationships  among  genes  in  a  given 
condition.  The  model  can  be  used  to  enhance  the 
understanding  of  gene  interactions  and  better  ways  of 
elucidating  environmental  and  drug-induced  effects. 

Large-scale  monitoring  of  gene  expression  such  as  DNA 
microarrays  [1-4]  is  considered  to  be  one  of  the  most 
promising  techniques  for  making  the  discovery  of  GRNs 
feasible  [5].  However,  the  task  of  inferring  GRNs  involves 
several  challenges  including  the  following:  (1)  the  number  of 
related  genes  is  very  large  compared  to  the  number  of  samples 
or  time  points,  (2)  the  observed  data  involve  a  significant 
amount  of  noise,  and  (3)  the  interaction  among  genes  displays 
complex  (nonlinear  and  dynamic)  relationships.  This 
challenge  provides  computer  scientists,  statisticians,  and 
engineers  with  opportunities  to  expand  their  knowledge  of 
intelligent  methods  to  provide  models  for  better  understanding 
of  biological  systems. 

The  field  of  system  modeling  plays  a  significant  role  in  the 
discovery  of  GRNs.  Several  system  modeling  approaches  have 
been  proposed  to  reverse-engineer  network  interactions 
including  a  variety  of  continuous  or  discrete,  static  or 
dynamic,  quantitative  or  qualitative  methods  [1-6]. 

The  use  of  computational  intelligence  (Cl)  methods  for 
system  modeling  has  gained  particular  interest,  because  they 
require  little  a  priori  knowledge  about  the  underlying  system 
and  the  model  can  be  derived  from  data.  Given  a  set  of  data, 


these  methods  discover  hidden  regularities  and  structures 
within  the  data.  Instead  of  requiring  patterns  to  be  known 
ahead  of  time,  they  search  automatically  for  patterns  that  are 
hidden  in  the  data. 

Several  Cl  methods  have  been  proposed  [7-14]  for  GRN 
reconstruction.  These  methods  have  been  successfully  applied 
in  both  cluster-  and  classification-based  approaches.  For 
example,  in  [9],  a  neural  model  is  used  to  simulate  the 
dynamics  of  the  lambda  phage  regulatory  system.  Middendorf 
et  al.  [15]  used  decision  trees  to  predict  whether  a  gene  is  up- 
or  down-regulated  in  a  particular  experiment  on  the  basis  of 
the  presence  of  binding  site  subsequences  (motifs)  in  the 
gene’s  regulatory  region  and  the  expression  level  of  regulators 
such  as  transcription  factors  in  the  experiment.  Soinov  et  al. 
[16]  applied  decision-tree  based  classifier  to  extract  simple 
rules  defining  gene  interrelations.  In  [17],  an  approach  is 
proposed  based  on  fuzzy  rules  of  a  known  activator/repressor 
model  of  gene  interaction.  This  algorithm  transforms 
expression  values  into  qualitative  descriptors  that  can  be 
evaluated  by  using  a  set  of  heuristic  rules  and  searches  for 
regulatory  triplets  consisting  of  activator,  repressor,  and  target 
gene.  This  approach,  though  logical,  is  a  brute  force  technique 
for  finding  gene  relationships.  It  involves  a  significant 
computation  time,  which  restricts  its  practical  usefulness. 
Also,  this  method  is  limited  to  the  study  of  the  interaction 
between  one  possible  positive  and  one  negative  regulator  for 
each  gene.  In  [18],  we  proposed  the  use  of  clustering  as  an 
interface  to  a  fuzzy  logic-based  method  to  improve  the 
computational  efficiency.  A  scalable  linear  variant  of  fuzzy 
logic  is  introduced  in  [19]  to  examine  the  interactions  of 
multiple  genes.  Genetic  algorithms  (GAs)  have  also  been 
applied  to  decipher  genetic  networks  from  gene  expression 
data  [10-12].  Shin  and  Iba  [12]  developed  an  inference 
algorithm  based  on  GAs  for  the  optimization  of  the  influence 
matrix  of  GRN.  In  [14],  GAs  and  ANNs  are  combined  to 
determine  gene  interactions  in  temporal  gene  expression  data. 

To  capture  the  nonlinear  and  dynamic  relationships,  we 
propose  to  model  GRNs  using  recurrent  neural  networks 
(RNNs),  which  consist  of  nonlinear  processing  elements 
(neurons)  that  possess  feedback  and  memory  units.  The 
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architecture  and  the  synaptic  weights  of  RNNs  are  optimized 
using  two  recently  introduced  swarm  intelligence  (SI) 
methods,  ant  colony  optimization  (ACO)  and  particle  swarm 
optimization  (PSO)  methods,  respectively.  The  hybrid  SI- 
RNN  algorithm  is  applied  to  infer  networks  of  interactions 
from  simulated  and  real-world  gene  expression  data,  yielding 
promising  results. 

The  paper  is  organized  as  follows.  Section  II  highlights  the 
steps  involved  in  using  Cl  for  system  modeling.  Section  III 
describes  our  proposed  SI-RNN  algorithm  for  inferring 
network  interactions.  Section  IV  presents  networks  inferred 
from  simulated  and  yeast  cell  cycle  gene  expression  data. 
Finally,  Section  V  concludes  the  paper. 

II.  Systems  Modeling  using  Cl  Methods 

Besides  selecting  a  suitable  Cl  paradigm,  developing  an 
intelligent  model  involves  four  steps:  data  preparation,  model 
structure  selection,  learning,  and  model  evaluation.  These 
steps  are  repeated  until  the  last  step  results  in  a  satisfactory 
performance. 

Many  times  the  “raw”  data  are  not  the  best  data  to  use  for 
modeling  a  Cl  paradigm.  Hence,  in  using  Cl  paradigms  to 
solve  real-world  problems,  it  is  important  to  transform  raw 
data  into  a  form  acceptable  to  the  paradigm.  The  first  step  is  to 
decide  what  the  inputs  and  outputs  are.  Inputs  that  are  not 
relevant  for  modeling  should  be  excluded.  The  next  step  is  to 
process  the  data  in  order  to  handle  missing  data,  remove 
outliers,  and  to  normalize  and  scale  the  data  into  acceptable 
range.  This  will  be  followed  by  model  structure  selection, 
which  includes  the  choice  of  neural  network  architecture, 
fuzzy  rules,  membership  functions,  fuzzy  operators,  genetic 
operators,  and  coding  scheme.  The  selection  of  neural  network 
architecture  includes  choosing  activation  functions, 
appropriate  number  of  layers,  number  of  neurons  in  each 
layer,  and  the  interconnection  of  the  neurons  and  the  layers. 
After  the  model  structure  is  selected,  its  free  parameters  will 
be  determined.  The  predominant  feature  of  Cl  paradigms  is 
that  they  learn  from  data.  We  define  learning  as  the  process  of 
finding  the  free  parameters  of  a  Cl  model  (e.g.  determining  the 
weights  of  a  neural  network).  After  learning  is  completed,  a  Cl 
paradigm  is  evaluated  for  its  performance  through  testing.  The 
purpose  of  this  testing  is  to  prove  the  adequacy  or  to  detect  the 
inadequacy  of  the  fitted  model.  The  latter  could  arise  from  an 
inappropriate  selection  of  network  topology,  too  small  or  too 
many  neurons,  or  from  insufficient  training  or  overtraining. 
Incorrect  input  node  assignments,  noisy  data,  error  in  the 
program  code,  or  several  other  effects  may  also  cause  a  poor 
fit.  The  aim  of  model  evaluation  is  to  insure  that  the  model  fit 
is  correct;  that  the  model  satisfies  the  desired  requirements, 
and  that  it  serves  as  a  general  model.  A  general  model  is  one 
whose  input-output  relationships  (derived  from  the  training 
dataset)  apply  equally  well  to  new  sets  of  data  (previously 
unseen  test  data)  from  the  same  problem  not  included  in  the 
training  set.  The  main  goal  of  intelligent  modeling  is  thus  the 
generalization  to  new  data  of  the  relationships  learned  on  the 


training  set. 

The  choice  of  appropriate  Cl  paradigm  is  critical  to  the 
modeling  process.  This  primarily  depends  on  the  complexity 
of  the  underlying  system  to  be  modeled,  the  available 
information  (a  priori  knowledge  and  data),  and  the  existence 
of  a  suitable  learning  algorithm.  In  recent  years,  ANNs  have 
been  employed  successfully  for  modeling  a  wide  range  of 
nonlinear  systems.  They  generally  consist  of  a  number  of 
interconnected  processing  elements  known  as  neurons.  The 
way  neurons  are  interconnected  or  how  the  inter-neuron 
connections  are  arranged  determines  the  architecture  of  a 
neural  network.  The  strengths  of  the  connections  (known  as 
weights  or  synaptic  weights)  are  adjusted  or  trained  to  achieve 
a  desired  overall  behavior  of  the  network.  The  most  popular 
architecture  is  a  feedforward  neural  network,  where  the 
neurons  are  grouped  into  layers.  All  connections  are 
feedforward;  that  is,  they  allow  information  transfer  only  from 
an  earlier  layer  to  the  next  consecutive  layers.  It  is  known  that 
ANN  can  sufficiently  approximate  the  nonlinear  mapping, 
learn  to  adapt  to  dynamics  of  uncertain  systems,  and  have 
strong  robustness  and  fault-tolerant  abilities  due  to  the  rich 
connection  and  nonlinear  activation  functions  of  the  neurons. 
In  light  of  the  above  advantages,  neural  network  based 
approaches  have  shown  the  superiority  over  the  well- 
established  and  proven  conventional  methods  for 
parameter/state  estimation  for  a  large  class  of  problems  and 
systems.  However,  in  order  to  perform  a  time  series  prediction 
or  build  a  model  of  a  dynamical  system  such  as  one  that 
represents  regulatory  interactions,  it  is  important  to  establish  a 
form  of  expansion  to  the  feedforward  neural  network,  so  that 
the  network  contains  some  type  of  memory  element.  This  can 
be  achieved  by  applying  time-delayed  inputs  to  feed  forward 
networks.  Alternatively,  an  RNN  can  be  built,  in  which  the 
outputs  of  some  neurons  are  fed  back  to  the  same  neurons  or 
to  other  neurons  in  the  network.  Thus,  signals  can  flow  in  both 
forward  and  backward  directions.  RNNs  have  a  dynamic 
memory  -  their  outputs  at  a  given  instant  reflect  the  current 
model  input  as  well  as  previous  inputs  and  outputs.  It  is  an 
ideal  candidate  for  modeling  dynamic  systems  such  as 
network  interactions  in  biological  systems. 

III.  Inferring  Network  Interaction  Using  SI-RNN 

In  building  an  RNN  to  infer  a  network  of  interactions,  the 
identification  of  the  correct  structure  and  determination  of  the 
free  parameters  (weights  and  biases)  to  mimic  measured  data 
is  a  challenging  task  given  the  limited  available  quantity  of 
data.  For  example,  in  inferring  a  GRN  from  microarray  data, 
the  number  of  time  points  is  considerably  low  compared  to  the 
number  of  genes  involved.  Considering  the  complexity  of  the 
biological  system,  it  is  difficult  to  adequately  describe  the 
pathways  involving  a  large  number  of  genes  with  few  time 
points.  In  this  paper,  we  apply  ACO  and  PSO  methods  to 
select  the  optimal  architecture  of  an  RNN  and  to  update  its 
free  parameters,  respectively.  ACO  is  a  discrete  optimization 
algorithm  that  has  been  successfully  used  for  combinatorial 


problems,  while  PSO  is  applied  for  continuous  optimization 
(i.e.,  the  variables  in  the  objective  function  can  assume  real 
values).  We  formulated  the  selection  of  RNN  structures  as  a 
combinatorial  problem  that  can  be  effectively  optimized  by 
ACO.  The  optimal  parameters  of  an  RNN  for  a  given  structure 
can  be  obtained  by  using  PSO  as  this  is  a  continuous 
optimization  problem. 

A.  Recurrent  Neural  Network 

Figure  la  shows  an  RNN,  where  the  output  of  each  neuron 
is  fed  back  to  its  input  after  a  unit  delay  and  is  connected  to 
other  neurons.  It  can  be  used  as  a  model  of  gene  regulatory 
network,  where  every  gene  in  the  network  is  considered  as  a 
neuron.  The  RNN  can  model  not  only  the  interactions  between 
genes  but  also  gene  self-regulation,  which  is  represented  by  a 
unit  delay  (z_1). 


Figure  1.  (a)  Architecture  of  a  fully  connected  RNN;  (b)  Details  of  a  single 
recurrent  neuron. 

Figure  lb  illustrates  the  details  of  the  ith  self- feedback 
neuron  (e.g.  ith  gene  in  the  GRN),  where  vh  known  as  the 
induced  local  field  (activation  level),  is  the  sum  of  the 
weighted  inputs  (the  regulation  of  other  genes)  to  the  neuron 
(z'th  gene);  and  cp(.)  represents  an  activation  function 
(integrated  regulation  of  the  whole  RNN  on  zth  gene),  which 
transforms  the  activation  level  of  a  neuron  into  an  output 
signal  (regulation  result).  The  induced  local  field  and  the 
output  of  the  neuron,  respectively,  are  given  by: 

vi(t)  =  fjwijxj(t-\)  +  bi  0) 

7=1 

x,(t)  =  (2) 

where  the  synaptic  weights  wn,  wi2,...,  define  the  strength 
of  connection  between  the  z'th  neuron  (e.g.  ith  gene)  and  its 
inputs  (e.g.  expression  level  of  genes).  Such  synaptic  weights 
exist  between  all  pairs  of  neurons  in  the  network.  bt  denotes 
the  bias  for  the  ith  neuron.  We  denote  s  as  a  structure  vector 
that  describes  the  architecture  of  the  network,  and  w  as  a 
weight  vector  that  consists  of  all  the  synaptic  weights  and 
biases  in  the  network,  s  and  w  are  adapted  during  learning 
to  yield  the  desired  network  outputs.  The  activation  function  cp 


introduces  nonlinearity  to  the  model.  When  information  about 
the  complexity  of  the  underlying  system  is  available,  a 
suitable  activation  function  can  be  chosen  (e.g.  linear,  logistic, 
sigmoid,  threshold,  hyperbolic  tangent  sigmoid  or  Gaussian 
function.)  If  no  prior  information  is  available,  our  algorithm 
uses  the  hyperbolic  tangent  sigmoid  function. 

As  a  cost  function,  we  use  the  mean-squared  error  between 
the  expected  output  and  the  network  output  across  time  (from 
the  initial  time  point  t0  to  the  final  time  point  tj)  and  across 
neurons  in  the  network.  The  cost  function  can  be  written  as: 

EW  =  TTr  0]2  (3) 

tfJS  t=tQ  i= 1 

where  xt(t)  and  x.(t)  are  the  true  and  predicted  values 
(expression  levels)  for  the  ith  neuron  (gene)  at  time  t.  The  goal 
is  to  determine  the  structure  vector  s  and  weight  vector  w 
that  minimize  this  cost  function.  We  propose  ACO  and  PSO  to 
optimize  s  and  w,  respectively. 

Note  that  the  above  RNN  is  self-evolutionary  and  can  be 
used  to  model  a  multi-step-ahead  prediction.  The  RNN  starts 
with  a  given  initial  condition,  evolves,  and  eventually  reaches 
final  states.  In  this  study,  we  inferred  a  network  from  a 
simulated  dataset  generated  by  a  five-node  network  using  both 
one-step-ahead  and  multi-step-ahead  forms,  which  resulted  in 
similar  performance.  For  real-world  gene  expression  data  that 
involve  noise  and  unequal  time  intervals,  we  use  the  one  step- 
ahead  prediction  method  during  training.  However,  the 
resulting  RNN  can  be  used  to  simulate  multiple-ahead 
predictions. 

B.  Ant  Colony  Optimization 

Ant  colony  optimization  studies  artificial  systems  that  take 
inspiration  from  the  behavior  of  real  ant  colonies.  The  basic 
idea  of  ACO  is  that  a  large  number  of  simple  artificial  agents 
are  able  to  build  good  solutions  to  solve  hard  combinatorial 
optimization  problems  via  low-level  based  communications. 

We  propose  to  use  ACO  to  optimize  the  structure  vector  s  . 
Each  possible  network  structure  s  is  defined  by  a 
combination  of  n  features  s  =  [sx  s2mmmSn],  where  Sj  is  an  n- 
bit  binary  string  that  indicates  which  neurons  are  controlled  by 
neuron  j.  Each  sj  is  selected  from  2n  candidate  features.  For 
each  neuron  j,  we  define  the  function  in  Eq.  (4)  to  determine 
the  probability  of  selecting  a  feature  i  among  the  T  candidate 
features: 

p/(k)=  ASIH  y  =  U«  /  =  2"  (4) 

yfr/w 

i=\ 

where  rj(k)  is  the  amount  of  pheromone  trail  for  the  ith 
feature  at  iteration  k.  At  k=  0,  r/  (/c)  is  set  to  a  constant  for  all 

features,  allowing  each  feature  to  have  equal  probability  of 
being  selected.  Thus,  in  the  first  iteration,  each  ant  chooses 
randomly  n  features  that  make  up  a  structure  (s  ,  a  trail).  Let 


s  be  an  ant  consisting  of  n  features  s  =[sx  s2m--sn]  . 
Depending  on  the  performance  of  s  ,  the  amount  of 
pheromone  trail  of  all  features  in  s  is  updated.  The 
performance  function  here  is  evaluated  on  the  basis  of 
mimicking  the  response  of  the  system  under  study.  To 
estimate  the  performance  of  s  ,  we  construct  an  RNN  that  has 
the  structure  defined  in  s  .  Then,  we  optimize  the  weights  of 
the  RNN  using  PSO.  The  response  of  the  resulting  RNN  will 
be  compared  with  the  measured/observed  response  of  the 
system  under  study.  The  amount  of  pheromone  trail  for  each 
element  in  s  is  updated  in  proportion  to  the  performance  of 
the  structure.  Assuming  that  the  /th  feature  for  the  yth  neuron 
was  in  s  ,  the  corresponding  amount  of  pheromone  trail  will 
be  updated  as  follows: 

zf(k  +  l)  =  p.z/(k)  +  Azl(k)  (5) 

where  p  is  a  constant  between  0  and  1,  representing  the 
evaporation  of  pheromone  trails,  and  A  T-(k)  is  an  amount 
proportional  to  the  performance  by  s  .  A  t-Qz)  is  set  to  zero,  if 
st€  s  .  This  update  is  made  for  all  N  ants  (sl9...,sN).  Note 

that  at  k=  0,  A i-(k)  is  set  zero  for  all  features.  The  updating  rule 
allows  trails  that  yield  good  performance  to  have  their  amount 
of  pheromone  trail  increased,  while  others  will  evaporate.  As 
the  algorithm  progresses,  features  with  large  amounts  of 
pheromone  trails  influence  the  probability  function  to  lead  the 
ants  towards  them. 

C.  Particle  Swarm  Optimization 

In  the  PSO  algorithm,  each  particle  is  represented  as  a 
vector  w.  and  instantaneous  trajectory  vector  Awi(k)  , 

describing  its  direction  of  motion  in  the  search  space  at 
iteration  k.  The  index  i  refers  to  the  /th  particle.  The  core  of 
the  PSO  algorithm  is  the  position  update  rule  (6)  which 
governs  the  movement  of  each  of  the  n  particles  through  the 
search  space. 

wL(k+X)  =  m>  (A) + Avt>  (A + 1) 

Ah;  (k+\)=  2(Ah;  (k)  +  O,  (wiheJk)  -  \\\  (/<:))+ <1>2  (w^Jk)  ~  (k))) 
where 
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At  any  instant,  each  particle  is  aware  of  its  individual  best 
position,  w.  best  (k) ,  as  well  as  the  best  position  of  the  entire 

swarm,  wG  hest  (k)  •  The  parameters  c\  and  c2  are  constants  that 

weight  particle  movement  in  the  direction  of  the  individual 
best  positions  and  global  best  positions,  respectively;  and  rXJ 
and  r2j- ,  y  =  l,2,...£>  are  random  scalars  distributed  uniformly 
between  0  and  1,  providing  the  main  stochastic  component  of 


the  PSO  algorithm. 

The  constriction  factor,  %  ,  may  also  help  to  ensure 
convergence  of  the  PSO  algorithm,  and  is  set  according  to  the 
weights  cx  and  c2  as  in  (7). 

2 

Z  =  I - ,  p  (p  =  cl+c2 ,  (p>  4  (7) 

\l-(p-yj(p2 -Aq>\ 

The  key  strength  of  the  PSO  algorithm  is  the  interaction 
among  particles.  The  second  term  in  (6),  0>2{wG^{k)-wi{k)) , 

is  considered  to  be  a  “social  influence”  term.  While  this  term 
tends  to  pull  the  particle  towards  the  globally  best  solution,  the 
first  term,  <&l  (w.  best  (k)  -  w?  (£)) ,  allows  each  particle  to  think  for 

itself.  The  net  combination  is  an  algorithm  with  excellent 
trade-off  between  total  swarm  convergence,  and  each 
particle’s  capability  for  global  exploration.  Moreover,  the 
relative  contribution  of  the  two  terms  is  weighted 
stochastically. 

The  algorithm  consists  of  repeated  application  of  the 
velocity  and  position  update  rules  presented  above. 
Termination  can  occur  by  specification  of  a  minimum  error 
criterion,  maximum  number  of  iterations,  or  alternately  when 
the  position  change  of  each  particle  is  sufficiently  small  as  to 
assume  that  each  particle  has  converged. 

Selection  of  appropriate  values  for  the  free  parameters  of 
PSO  plays  an  important  role  in  the  algorithm’s  performance. 
In  our  study,  parameters  cx  and  c2  were  arbitrarily  selected  (cx 
=  2.05,  c2  =  2.05),  with  constriction  factor,  /,  determined  by  (7) 
and  the  maximum  velocity  was  set  at  2. 


D.  SI-RNN 

In  this  section,  we  illustrate  how  the  two  SI  methods,  ACO 
and  PSO,  work  together  to  optimize  both  s  and  w  of  an 
RNN  to  mimic  the  response  of  an  unknown  network  of 
interactions.  Each  node  in  the  true  network  is  represented  by  a 
neuron  in  the  RNN.  We  assume  that  the  number  of  nodes  in 
the  network  is  known  (e.g.  the  number  of  genes  for  which  a 
gene  regulatory  network  is  to  be  modeled  is  known),  but  the 
way  the  nodes  interact  is  assumed  to  be  unknown. 

s)  j  s2 

J  =[0  11001000] 

1  ^  l 

Vp  =  [  0.7  0.8  0.15] 


Figure  2.  True  network  (a);  randomly  selected  network  (b). 


ACO  starts  with  initial  candidate  ants  that  define  various 
structures.  For  example,  let  Fig.  2(a)  be  the  structure  for  the 
true  interaction  that  has  three  nodes  and  Fig.  2(b)  be  one  of  the 
randomly  selected  initial  structures  by  ACO.  For  each  of  these 


two  systems,  the  corresponding  structure  vector  s  is  shown  in 
the  figure,  where  Sj  (j  =  1,  2,  3)  is  a  three-bit  binary  string  that 
indicates  which  neurons  (including  itself)  are  controlled  by  the 
zth  neuron.  For  example,  if  s\  =  [0  1  1],  it  implies  that  the  first 
neuron  controls  all  others  except  itself  and  s2  =  [0  0  1]  implies 
that  the  second  neuron  controls  the  third  neuron  only. 

PSO  searches  for  the  optimal  weight  vector  w  to  minimize 
the  difference  between  the  output  of  the  true  network  and  the 
RNN  using  the  training  data.  Only  the  elements  of  w  that 
correspond  to  nonzero  entries  in  s  are  updated  by  PSO.  For 
example,  the  weight  vector  w  in  Fig.  2b  contains  only  two 
variables:  0.25  and  0.4,  corresponding  to  the  two  nonzero 
entries  in  the  structure  vector  s  . 

The  optimal  weight  vectors  for  all  randomly  selected 
structures  (ants)  are  tested  with  previously  unseen  validation 
data.  The  performance  of  each  particle  in  simulating  the 
validation  data  is  returned  to  ACO  to  update  the  trails  of  the 
ants  in  the  search  space  (i.e.,  update  the  structure  vector  s  ). 
The  new  s/s  in  s  are  used  to  construct  a  new  candidate 
structure.  This  will  lead  to  a  structure  that  is  more  similar  to 
the  global  best  structure  than  the  previous  one.  The 
assumption  in  this  algorithm  is  that  the  prediction  error  of  an 
arbitrary  network  will  be  larger  than  a  network  that  matches 
the  correct  structure.  Through  subsequent  iterations,  the  ACO 
and  PSO  search  for  the  optimal  structure  vector  s  and  weight 
vector  w  to  make  accurate  predictions. 

|  Raw  Data 
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Figure  3.  Flowchart  for  the  SI-RNN  algorithm. 

Figure  3  illustrates  the  overall  algorithm,  which  involves 
the  two  major  components:  (i)  ACO  that  initially  generates 
various  topologies  randomly  and  selects  the  most  optimal 


structure  iteratively  based  on  the  prediction  performance  of 
the  topologies  in  a  validation  dataset  and  (ii)  PSO  that 
determines  the  free  parameters  of  a  given  topology  with  the 
aim  of  minimizing  the  error  between  the  training  data  and  the 
outputs  of  the  networks  whose  topologies  are  determined  by 
ACO.  The  performance  of  the  PSO  particles  in  predicting  the 
validation  dataset  is  used  to  iteratively  update  the  network 
structure.  The  process  continues  until  satisfactory  result  is 
found  or  maximum  number  of  iterations  is  reached.  The  final 
network  will  be  evaluated  via  independent  testing  dataset. 

IV.  Data  and  Results 

The  SI-RNN  algorithm  is  evaluated  in  inferring  networks  of 
interactions  from  artificial  and  yeast  cell  cycle  data. 

A.  Artificial  Data 

We  applied  the  SI-RNN  approach  to  identify  the  network  in 
Fig.  4.  We  generated  three  datasets  (training,  validation,  and 
testing)  with  different  initial  conditions.  Each  dataset 
consisted  of  20  time  points.  These  artificial  data  sets  are 
created  to  ascertain  the  ability  of  the  algorithm  to  “rediscover” 
the  underlying  network  that  generated  the  data. 

An  RNN  model  of  five  neurons  with  hypothetic  tangent 
sigmoid  activation  function  was  trained  using  the  SI-RNN 
algorithm.  PSO  used  the  training  dataset  to  determine  the 
optimal  weight  vector  w  for  each  structure  vector  s  defined 
by  ACO.  The  performance  of  each  structure  in  predicting  the 
outputs  of  the  network  in  the  validation  dataset  is  used  by 
ACO  to  determine  the  optimal  structure.  The  algorithm  was 
run  100  times.  In  each  run,  Eq.  3  was  evaluated  1000  times  to 
identify  the  structure  that  leads  to  the  least  cost.  54  runs  (out 
of  100)  predicted  a  RNN  with  identical  structure  to  Fig.  4. 
16%  of  the  runs  also  predicted  the  true  network  structure,  but 
they  had  two  or  three  additional  connections.  The  remaining 
30%  consisted  of  arbitrary  structures.  Fig.  5  shows  the  outputs 
of  the  true  network  and  the  predicted  RNN  for  the  testing 
dataset. 

3q(^)=0.5v5(^-l) 
x2(t)=03x[(t-T) 
x3(t)  =-0.6x2(/-1)  +x4(t  -1) 
xA{t)=0.l[xA{t  -l))2 
x5(0^.8^M 

Figure  4.  A  simulated  five-node  network. 

B.  Yeast  Cell  Cycle  Data 

The  yeast  cell  cycle  data  collected  by  Spellman  et  al.  [20] 
consist  of  six  time  series  (cln3,  clb2,  alpha,  cdcl5,  cdc28  and 
elu)  expression  measurements  of  the  transcript  (mRNA)  levels 
of  S.  cerevisiae  genes. 

To  test  our  approach,  we  chose  five  cyclin  genes  (CLB1, 
CLB2,  CLB5,  CLB  6  and  CDC28),  which  are  involved  in  cell- 
cycle  regulation.  The  expression  levels  of  these  genes  in  three 
time  series  measurements  were  used  to  construct  a  GRN.  In 


these  datasets,  the  biological  samples  were  synchronized  by 
three  different  methods:  a  factor  arrest,  arrest  of  a  cdcl5,  and 
cdc28  temperature-sensitive  mutant.  We  used  the  cdcl5 
dataset,  which  has  24  experimental  conditions,  as  training 
dataset.  The  other  two  datasets,  alpha  and  cdc  28,  which  have 
18  and  17  time  points,  were  used  as  validation  and  testing 
datasets,  respectively. 
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Figure  5.  Original  and  predicted  outputs  of  the  testing  dataset. 
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Figure  6.  Original  and  predicted  outputs  of  five  genes  in  the  testing  dataset 
(cdc28  dataset). 


PSO  used  the  training  dataset  to  determine  the  optimal 
weight  vector  w  for  each  structure  vector  s  defined  by 
ACO.  The  performance  of  each  structure  in  predicting  the 
outputs  of  the  network  in  the  validation  dataset  is  used  by 
ACO  to  determine  the  optimal  structure.  The  algorithm  was 
run  10  times.  In  each  run,  Eq.  3  was  evaluated  1000  times  to 
identify  the  structure  that  leads  to  the  least  MSE.  To  improve 
the  prediction  accuracy,  only  the  connections  obtained  in  at 
least  50%  of  the  runs  were  selected.  Fig.  6  shows  the 
measured  gene  expression  data  and  the  outputs  of  the 
predicted  RNN  for  each  gene  in  the  testing  dataset. 

As  the  true  GRN  that  governs  the  interaction  among  the  five 
genes  is  not  available,  the  accuracy  of  the  network  is 
determined  by  how  well  it  fits  the  measured  gene  expression 
data.  To  get  insight  into  the  performance  of  the  SI-RNN 
method,  we  randomly  generated  100  RNN  structures  and 
optimized  their  parameters  using  PSO.  The  average  and 
standard  deviation  of  the  MSE  for  the  100  randomly  generated 
RNN  were  0.25  and  0.14,  respectively,  while  the  optimal 
RNN  structure  found  through  the  SI-RNN  method  yielded  an 
average  MSE  =  0.1  and  standard  deviation  =  0.08  in  10  runs. 
Note  that  the  MSE  is  calculated  for  the  data  normalized 
between  0  and  1. 

To  further  evaluate  the  GRN  identified  by  SI-RNN,  we 
searched  for  known  interactions  in  a  S.  cerevisiae  database 
provided  by  the  PathwayStudio  software  (Ariadne  Genomics, 
Rockville,  MD).  Figures  7a  and  7b  depict  the  interactions 
among  the  five  genes  found  by  the  PathwayStudio  software 
and  by  our  hybrid  SI-RNN  algorithm,  respectively.  Table  I 
illustrates  the  gene  relationships  depicted  in  Fig.  7. 

Among  the  1 1  relations  predicted  by  SI-RNN,  five  of  them 
concur  with  the  known  interactions  obtained  by 
PathwayStudio.  Two  gene  self-regulations  (CLB1  and  CLB2) 
were  found  by  SI-RNN,  but  not  by  PathwayStudio.  Four  other 
relations  (presented  by  dotted  lines  in  Fig.  7b)  are  also  found 
by  PathwayStudio,  but  with  reversed  direction  of  regulation. 
PathwayStudio  searches  for  known  interactions  on  the  basis  of 
literature  mining  through  natural  language  processing 
methods.  From  the  literature  used  by  PathwayStudio,  it 
appears  that  the  three  “reversed”  relations  are  described  to 
have  some  relations  without  a  defined  direction  of  regulation. 
For  example,  CFB5  and  CFB6  are  essential  for  sporulation 
because  they  are  required  for  premeiotic  DNA  replication 
[21].  This  indicates  that  there  maybe  some  relationship 
between  these  two  genes,  but  not  necessarily  a  directed 
regulation  from  CFB5  to  CFB6.  Another  example  is  the 
relationship  between  CFB6  and  CDC28.  In  [22],  it  is  stated 
that  the  actual  initiation  events  require  the  activities  of  at  least 
two  protein  kinases,  the  cyclin-dependent  kinase  (CDK) 
Cdc28p  associated  with  cyclin  B  (Clb5p  or  Clb6p)  and  the 
Cdc7p  kinase  associated  with  its  regulatory  subunit  Dbf4p.  No 
direct  regulation  information  is  provided. 


TABLE  I 

Known  Relations  among  Five  Genes  from  PathwayStudio  Software. 


Relation  Type 

Symbol 

Predicted  by  SI-RNN 

Expression 

CLB1  <—  CLB6 

yes  (reversed) 

Expression 

CLB1  <—  cdc28 

no 

Expression 

CLB2  <—  cdc28 

no 

Expression 

CLB1  <+-  CLB2 

yes 

Regulation 

CLB6  — >  CLB5 

yes  (reversed) 

Regulation 

CLB6  -+>  CLB2 

no 

Regulation 

CLB1  <+-  CLB5 

yes  (reversed) 

MolSynthesis 

CLB1  -+>  CLB2 

yes 

Genetic  Interaction 

CLB2  — -  cdc28 

no 

Direct  Regulation 

CLB6  — +>  cdc28 

yes  (reversed) 

Direct  Regulation 

CLB5  -+>  cdc28 

yes 

Direct  Regulation 

CLB2  -+>  cdc28 

yes 

Direct  Regulation 

CLB2  <+-  CLB5 

no 

Direct  Regulation 

CLB1  -+>  cdc28 

yes 

(a) 

I - 1 


Figure  7.  Cell  cycle  pathway  in  five  genes:  CDC28,  CLB1,  CLB2,  CLB5  and 
CLB6.  (a)  Known  interactions  found  by  the  PathwayStudio  software,  (b) 
Result  of  the  proposed  hybrid  SI-RNN  method:  correctly  identified  relations 
use  the  same  line  as  (a),  reversed  relations  use  dotted  lines,  and  two  additional 
self-regulations  are  indicated  (CLB1  and  CLB2). 


V.  Conclusion 

In  this  paper,  we  explored  the  combined  advantages  of  the 
nonlinear  and  dynamic  properties  of  RNN,  and  the  global 
search  capabilities  of  swarm  intelligence  methods  to  infer 
network  interactions.  We  evaluated  the  performance  of  the 
algorithm  using  data  generated  from  a  simulated  five-node 
network,  and  time  course  gene  expression  data  for  five  yeast 
genes.  Although  the  algorithm  yielded  promising  result  in 
predicting  the  simulated  network,  due  to  the  stochastic 
properties  of  the  algorithms,  not  all  runs  identify  the  correct 
network  structure.  Our  future  work  will  focus  on  improving 
the  rate  at  which  the  correct  structure  is  identified. 

In  inferring  the  GRN  that  govern  the  five  genes,  about  half 
of  the  relations  predicted  by  SI-RNN  were  verified  by 
published  literature.  Three  additional  interactions  are  also 
verified  in  the  literature  to  have  some  relations,  although  no 
specific  direction  of  relation  was  given.  However,  these 
relations  predicted  by  SI-RNN  were  reversed  compared  to  the 
output  of  the  PathwayStudio  software.  Two  self-regulations 
were  predicted  by  SI-RNN,  but  they  could  not  be  verified  by 
the  PathwayStudio  software. 

Due  to  noise  and  insufficient  time  points  in  real  gene 
expression  data,  we  anticipate  challenges  in  applying  the 
proposed  hybrid  SI-RNN  method  to  infer  gene  regulatory 
networks  that  involve  large  number  of  genes.  To  address  these 
challenges,  we  plan  to  incorporate  known  gene  interactions 
and  genomic  information  into  the  SI-RNN  algorithm. 
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Summary 

A  serious  obstacle  to  successful  treatment  of  estrogen  receptor  (ER)-positive  human  breast  cancer  is  cell  resistance  to 
tamoxifen  (TAM)  therapy.  Here  we  show  that  the  electrophile  disulfide  benzamide  (DIBA),  an  ER  zinc  finger  inhibitor,  blocks 
ligand-dependent  and  -independent  cell  growth  of  TAM-resistant  breast  cancer  in  vitro  and  in  vivo.  Such  inhibition  depends 
on  targeting  disruption  of  the  ER  DNA-binding  domain  and  its  communication  with  neighboring  functional  domains, 
facilitating  ERa  dissociation  from  its  coactivator  AIB1  and  concomitant  association  with  its  corepressor  NCoR  bound  to 
chromatin.  DIBA  does  not  affect  phosphorylation  of  HER2,  MAPK,  AKT,  and  AIB1,  suggesting  that  DIBA-modified  ERa 
may  induce  a  switch  from  agonistic  to  antagonistic  effects  of  TAM  on  resistant  breast  cancer  cells. 


Introduction 

The  selective  estrogen  receptor  modulator  (SERM)  tamoxifen 
(TAM),  which  binds  to  the  estrogen  receptor  a  (ERa)  and  partially 
inhibits  its  activity,  is  the  most  prolific  therapeutic  drug  for  the 
treatment  of  ER-positive  breast  cancer  (Osborne,  1998).  Adju¬ 
vant  therapy  studies  of  TAM  show  a  40%-50%  reduction  in 
the  odds  of  recurrence  and  reduced  mortality.  Unfortunately, 
advanced  breast  cancers  that  initially  respond  well  to  TAM 
eventually  become  refractory  to  this  compound  (McDonnell 
and  Norris,  2002;  Jordan,  2004;  Osborne  et  al.,  2003;  Shou 
et  al.,  2004). 

ER  functions  in  the  nucleus  as  a  transcriptional  regulator  of 
specific  genes  (Tsai  and  O’Malley,  1994).  The  structural  organi¬ 
zation  of  ERa  consists  of  a  ligand-independent  transcription- 
activation  domain  (AF-1  domain),  a  DNA-binding  domain  (DBD), 
a  ligand-binding  domain  (LBD),  and  a  ligand-dependent  trans¬ 
activation  domain  (AF-2  domain)  (Kumar  et  al.,  1987;  Ruff 


et  al.,  2000).  Estrogen  binding  to  ER  alters  its  conformation,  trig¬ 
gers  receptor  dimerization,  and  directly  facilitates  binding  of  the 
receptor  complex  to  promoter  regions  of  target  genes,  including 
sites  known  as  estrogen-responsive  elements  (ERE),  or  indi¬ 
rectly  through  transcription  factors  such  as  AP-1  (Kushner 
et  al.,  2000).  The  recruitment  of  coactivators  such  as  AIB1  and 
other  proteins  with  acetyltransferase  activity  helps  to  unwind 
the  chromatin,  allowing  transcription  to  occur  (Brzozowski 
et  al.,  1997;  Glass  and  Rosenfeld,  2000;  Shang  et  al.,  2000; 
Shiau  et  al.,  1998;  Smith  et  al.,  1997).  In  contrast,  the  ER  con¬ 
formation  induced  by  the  binding  of  SERMs  like  TAM  favors 
the  recruitment  of  corepressors  NCoR/SMRT  and  deacetylases 
that  inhibit  transcriptional  activity  in  TAM-sensitive  breast  can¬ 
cer  cells  (Keeton  and  Brown,  2005;  Kurebayashi  et  al.,  2000; 
Mak  et  al.,  1999;  Osborne  et  al.,  2003;  Shou  et  al.,  2004).  How¬ 
ever,  acquired  resistance  can  be  caused  by  alterations  in  the  ER 
signal  transduction  pathway,  converting  the  inhibitory  SERM- 
ERa  complex  to  a  growth  stimulatory  signal  (Jordan,  2004). 


SIGNIFICANCE 

Acquired  resistance  to  antiestrogens  is  a  major  challenge  to  the  clinical  management  of  initially  endocrine-responsive  metastatic 
breast  cancer.  We  have  previously  found  that  electrophilic  DIBA  and  benzisothiazolone  derivatives  inhibited  TAM-sensitive  breast  can¬ 
cer  cells  by  preferentially  disrupting  the  vulnerable  zinc  fingers  within  the  ER  DNA-binding  domain.  Here  we  describe  how  DIBA  restores 
the  antagonistic  action  of  TAM  in  resistant  breast  cancer  cells  through  targeted  disruption  of  the  ER  DNA-binding  domain  and  its  inter¬ 
action  with  the  proximal  N-terminal  domain  to  suppress  ligand-dependent  and  -independent  ER  transcription  and  influence  the  recruit¬ 
ment  of  cofactor  to  the  ER.  These  results  show  that  small-molecule  modification  of  the  ER  zinc  finger  may  alter  coactivator/corepressor 
functions,  which  are  particularly  relevant  to  TAM  resistance. 


CANCER  CELL  10,  487-499,  DECEMBER  2006  ©2006  ELSEVIER  INC.  DOI  1 0.1 01 6/j.ccr.2006.09.01 5 


487 


ARTICLE 


Growing  evidence  indicates  that  crosstalk  between  ER  and 
growth  factor  receptor  signaling  pathways  (Brockdorff  et  al., 
2003;  Ibrahim  and  Yee,  2005;  Osborne  et  al.,  2005),  especially 
the  insulin-like  growth  factor  receptor  (IGFR)  family  and  the  epi¬ 
dermal  growth  factor  receptor  (EGFR)  family  (such  as  cErbB2 
[FIER2]),  is  one  of  the  mechanisms  for  resistance  to  endocrine 
therapy  in  breast  cancer  (Schiff  et  al.,  2004).  In  tumors  with 
abundant  ER,  AIB1 ,  and  HER2,  TAM  behaves  as  an  ER  agonist 
and  stimulates  tumor  growth  (Osborne  et  al.,  2005).  High  levels 
of  activated  AIB1  could  reduce  the  antagonist  effects  of  TAM, 
especially  in  tumors  that  also  overexpress  the  HER2  receptor 
that  activates  MAPKs.  TAM  resistance  may  also  be  produced 
by  decreased  levels  of  the  corepressor  NCoR  (Fujita  et  al., 
2003;  Lavinsky  et  al.,  1998;  Osborne,  1998). 

The  ER-DBD  contains  two  nonequivalent  Cys4  zinc  fingers 
(Laity  et  al.,  2001 ;  Ruff  et  al.,  2000;  Schoenmakers  et  al.,  1999; 
Wikstrom  et  al.,  1 999),  which  function  cooperatively  in  ER  dimer¬ 
ization  and  DNA  binding  by  stabilizing  the  secondary  and  tertiary 
structure  of  the  ER-DNA  complex  (Maynard  and  Coveil,  2001; 
Predki  and  Sarkar,  1992;  Schwabe  et  al.,  1993),  leading  to 
ligand-dependent  ER  transactivation  and  ER-mediated  breast 
cancer  cell  and  tumor  growth.  Moreover,  interdomain  commu¬ 
nication  between  the  N-terminal  AF-1  domain  and  DBD  of 
the  nuclear  receptors  helps  modulate  structure-  and  ligand- 
independent  functions  of  receptors  (Brodie  and  McEwan, 
2005;  Kumar  and  Thompson,  2003;  Shao  et  al.,  1 998;  Takimoto 
et  al.,  2003).  We  have  previously  found  that  electrophilic  DIBA 
and  benzisothiazolone  derivatives  produced  anticancer  activity 
in  TAM-sensitive  human  breast  cancer  cells  by  preferentially 
disrupting  the  vulnerable  ER  zinc  fingers,  thus  blocking  ER 
DNA  binding  and  transactivation  (Wang  et  al.,  2004).  Since  this 
anti-breast-cancer  strategy  targeted  ER  at  the  level  of  its  DNA 
binding,  rather  than  the  classical  antagonism  of  estrogen  bind¬ 
ing,  it  is  relevant  to  explore  whether  DIBA  has  the  capacity  to 
inhibit  the  growth  of  TAM-resistant  breast  cancer  cells. 

In  this  report,  we  investigated  how  DIBA  restored  the  antago¬ 
nist  action  of  TAM  on  breast  cancer,  which  was  dependent  on 
targeting  disruption  of  the  ER  DNA-binding  domain  and  its  com¬ 
munication  with  neighboring  transcription  domains.  Moreover, 
DIBA  reduced  ER  association  with  coactivator  AIB1  and  en¬ 
hanced  ER  association  with  corepressor  NCoR.  These  findings 
provided  the  proof  of  principle  for  a  potential  for  DIBA  applicable 
to  TAM-resistant  breast  cancer. 


Results 

DIBA  suppresses  TAM-resistant  breast  cancer 
cell  growth 

First  we  explored  whether  DIBA  affects  estrogen-mediated 
growth  of  TAM-resistant  breast  cancer  cells.  MCF-7/LCC2  is 
a  selective  ER-positive,  TAM-resistant  cell  line  (Brunner  et  al., 
1993;  Lilling  et  al.,  2000).  The  specific  ER  ligand  17|3-estradiol 
(E2)  stimulated  [3H]thymidine  incorporation  in  MCF-7/LCC2 
and  its  parent  MCF-7  cells,  but  the  degree  of  stimulation  in 
MCF-7/LCC2  is  significantly  less  than  that  observed  in  E2- 
treated  MCF-7  cells  (Figures  1A  and  IB).  4-Hydroxytamoxifen 
(4-OH-TAM)  significantly  inhibited  MCF-7  cells,  with  an  ED50 
of  0.1  jiM.  A  low  dosage  of  DIBA  enhanced  TAM  sensitivity, 
the  ED50  deceasing  2-fold  (0.05  jiM)  (Figure  1A).  The  TAM-resis¬ 
tant  cell  line  MCF-7/LCC2  validated  with  relative  resistance; 


however,  a  small  dosage  of  DIBA  (5  jiM)  restored  4-OH-TAM 
sensitivity,  achieving  over  90%  inhibition  of  E2-driven  prolifera¬ 
tion  at  the  lowest  dosage  tested  of  4-OH-TAM  (0.05  jiM). 
Similarly,  DIBA  inhibited  cell  proliferation  of  MCF-7/HER2-18 
(Figure  1C),  another  TAM-resistant  MCF-7  derivative  engi¬ 
neered  to  overexpress  HER2  (Benz  et  al.,  1993),  and  different 
types  of  ER-positive  and  TAM-resistant  breast  carcinoma  cell 
lines  including  BT474  (Figure  ID),  which  expresses  ER  and  is 
naturally  gene  amplified  for  HER2  and  AIB1  (Lin  et  al.,  1990; 
Anzick  et  al.,  1997),  and  epithelial  ZR-75  cells  (Figure  IE) 
(Hoffmann  et  al.,  2004)  in  a  dose-dependent  manner.  These 
observations  suggested  that  DIBA  effectively  restored  the 
antagonist  action  of  TAM  on  growth  of  TAM-resistant  breast 
cancer  cells. 

In  TAM-resistant  cells,  peptide  growth  factor  signaling  path¬ 
ways  appear  to  be  important  in  modifying  cell  behavior,  growth, 
and  survival  (Brockdorff  et  al.,  2003;  Ibrahim  and  Yee,  2005). 
Therefore,  we  examined  whether  DIBA  impacted  TAM-resistant 
cell  growth  mediated  by  stimulation  of  exogenous  peptide 
growth  factors.  MCF-7/LCC2  cells  (Figure  IF)  were  stimulated 
by  IGF-1  alone  or  IGF-1  plus  4-OH-TAM.  TAM  did  not  block 
IGF-1 -driven  cell  proliferation.  However,  adding  DIBA  at  even 
1  jiM  was  sufficient  to  restore  TAM  inhibitory  functions.  These 
data  demonstrated  that  DIBA  also  suppressed  TAM-resistant 
cell  growth  mediated  by  growth  factors. 

Efficacy  of  DIBA  on  TAM-resistant  breast  cancer 
tumor  growth  in  vivo 

The  in  vivo  efficacy  of  the  DIBA  was  tested  using  nude  mice 
bearing  human  MCF-7/LCC2  breast  carcinoma  xenografts. 
4-OH-TAM  alone  did  not  significantly  affect  tumor  growth. 
DIBA  alone  resulted  in  a  dose-dependent  inhibition  of  tumor 
growth,  and  a  high  dose  (30  mg/kg)  of  DIBA  reduced  tumor  vol¬ 
ume  to  almost  50%.  Moreover,  treatment  with  4-OH-TAM  plus 
DIBA  diminished  tumor  to  undetectable  levels  (Figure  2A).  Histo- 
pathological  analysis  (Figure  2B)  showed  a  typical  hypercellular 
solid  carcinoma  invading  the  dermis  and  subcutaneum,  and  the 
tumor  cells  had  a  high  nuclear  grade  with  frequent  mitosis  in  the 
control  vehicle  (upper  panel)  or  4-OH-TAM  alone-treated  mice 
(middle  panel).  In  contrast,  marked  reduction  in  tumor  volume, 
partial  encapsulation  by  fibrous  connective  tissue,  and  no  signif¬ 
icant  invasion  into  surrounding  skin  tissue  were  observed  in  the 
mice  treated  with  DIBA  plus  4-OH-TAM  (lower  panel).  These 
tumor  cells  with  a  low  nuclear  grade,  focal  glandular  differentia¬ 
tion,  and  no  frequent  mitosis  or  necrosis  were  seen  under  higher 
magnification.  No  apparent  toxicity  was  observed  in  liver  or 
kidney  in  Dl BA-treated  mice,  nor  were  there  any  significant 
changes  in  body  weight  gain  compared  with  control  mice 
(data  not  shown).  Therefore,  the  data  demonstrate  that  DIBA 
effectively  reduces  the  growth  of  MCF-7/LCC2  TAM-resistant 
tumors  in  mice. 

Synergism  between  DIBA  and  TAM 
on  cell-cycle  progression 

Using  propidium  iodide  (PI)  staining  and  fluorescence-activated 
cell  sorting  (FACS)  analysis,  we  further  evaluated  TAM-treated 
cells  within  the  cell  cycle  in  the  presence  of  DIBA  (Figure  2C). 
E2-treated  MCF-7/LCC2  cells  showed  decreased  cells  in  the 
G0/G1  phase  and  an  increased  percentage  of  cells  in  the  S 
and  G2/M  phases.  Cells  treated  with  TAM  had  a  weak  inhibitory 
effect  on  E2,  increasing  the  percentage  of  cells  in  S/G2/M.  By 
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Figure  1.  DIBA  is  a  potent  inhibitor  of  TAM-resistant  breast  cancer  cell  proliferation 

A-E:  Proliferation  of  MCF-7  (A),  MCF-7/LCC2  (B),  MCF-7/HER2-1 8  (C),  BT474  (D),  or  ZR-75  (E)  cells  was  examined  by  [3H]thymidine  incorporation  assay.  Starved 
cells  were  treated  with  DIBA  for  2  hr,  stimulated  with  (filled  bars)  or  without  (hatched  bars)  50  nM  E2,  incubated  with  increasing  concentrations  of  4-OH-TAM, 
and  analyzed  48  hr  later.  Data  shown  represent  mean  ±  SEM. 

F:  Proliferation  of  starved  MCF-7/LCC2  cells  induced  by  50  nM  4-OH-TAM  or  5  ng/ml  IGF- 1  was  also  examined  after  treatment  with  increasing  concentrations  of 
DIBA.  Data  shown  represent  mean  ±  SEM. 


contrast,  in  cells  cotreated  with  E2  plus  IC1 1 82780,  the  changes 
in  cell-cycle  status  and  growth  induced  by  E2  were  significantly 
inhibited.  In  the  presence  of  DIBA  combined  with  TAM,  cell- 
cycle  phase  distribution  induced  by  E2  shows  a  significant  in¬ 
crease  (from  60.7%  to  81 .7%)  of  cells  in  the  G0/G1  phase,  a  de¬ 
crease  (from  30.1  %  to  6.4%)  in  the  S  phase,  a  decrease  (from 
7.9%  to  5.4%)  in  the  G2/M  phase,  and  an  increase  (from  1 .2% 
to  6.5%)  in  the  sub-GI  phase.  Also,  DIBA  enhanced  the  inhibi¬ 
tory  effect  of  ICI  182780  on  E2-stimulated  cell  growth.  The 
FACS  data  further  confirmed  that  DIBA  restored  the  antagonist 
action  of  TAM  on  cell  proliferation  of  TAM-resistant  breast  can¬ 
cer  cells  analyzed  by  the  [3H]thymidine  incorporation  assay. 


ER  is  necessary  for  synergism  between  DIBA  and  TAM 

To  determine  whether  targeted  disruption  of  ER  is  necessary  for 
DIBA  to  suppress  cell  growth  of  TAM-resistant  breast  cancer 
cells,  we  used  BT474,  an  ER-positive  but  TAM-resistant  breast 
cancer  cell  line,  as  a  model  system  to  examine  the  effect  of  de¬ 
pletion  of  ER  on  DIBA  inhibition  of  cell  growth  of  TAM-resistant 
cells  (Figure  3A).  The  ER  expressed  in  these  cells  was  knocked 
down  by  using  ERa-siRNA.  The  decreased  level  of  ER  was 
confirmed  by  western  blot  (Figure  3A,  inset).  Under  identical 
conditions,  DIBA  rendered  TAM  inhibition  on  parent  ER-positive 
cells,  but  was  not  able  to  sensitize  TAM’s  suppression  of  growth 
of  ER-depleted  breast  cancer  cells.  These  data  suggest  that 
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Figure  2.  Synergism  inhibition  between  DIBA  and 
TAM  on  in  vivo  tumor  growth  and  cell-cycle  pro¬ 
gression 

A:  Dose-dependent  effect  of  DIBA  and  4-OH- 
TAM  on  growth  of  MCF-7/LCC2  tumor  in  mice. 
Data  shown  represent  mean  ±  SEM  (n  =  10 
mice  per  group). 

B:  Morphology  of  MCF-7/LCC2  tumors  treated 
with  vehicle  (upper  panel),  4-OH-TAM  at 
20  mg/kg/day  (middle  panel),  or  4-OH-TAM  at 
20  mg/kg/day  plus  DIBA  at  30  mg/kg  (lower 
panel). 

C:  MCF-7/LCC2  cells  were  synchronized  by 
serum  starvation,  pretreated  with  5  [xM  DIBA, 
and  then  stimulated  with  50  nM  E2,  50  nM  TAM, 
or  1  |xM  ICI  182780.  Cell-cycle  distribution  was 
examined  by  PI  staining  and  FACS  analysis.  The 
results  represent  three  independent  experiments 
(mean  ±  SEM). 


inhibition  of  DIBA  on  growth  of  TAM -resistant  breast  cancer 
cells  depends  on  ER. 

DIBA  inhibits  ER  binding  to  DNA 

To  clarify  whether  the  DIBA  alters  estrogen-  orTAM-bound  ER’s 
ability  to  bind  to  its  cognate  ERE  in  TAM-resistant  breast  cancer 
cells,  we  performed  electrophoretic  mobility  shift  assays  (EMSA) 
using  nuclear  extracts  obtained  from  MCF-7/LCC2  cells  (Fig¬ 
ure  3B).  E2-  or  4-OH-TAM-treated  cells  displayed  considerable 
ERE  DNA-binding  complexes,  which  could  be  partially  super- 
shifted  with  anti-ER,  but  not  normal  rabbit  serum,  confirming 
the  specificity  of  these  binding  complexes.  DIBA  significantly 
decreased  (80%)  the  E2-  or  TAM-induced  ERE  DNA-binding 
activity.  In  a  similar  experiment  on  androgen  receptor  (AR)  in 
MCF-7/LCC2  (Figure  3C),  DIBA  did  not  inhibit  AR  DNA-binding 
activity. 

Next,  we  examined  whether  DIBA  affects  ER  binding  to 
probes  containing  the  AP-1 ,  a  nontypical  repeat  element.  Estro¬ 
gen  or  TAM  induced  substantial  AP-1  binding  activity  (Fig¬ 
ure  3D).  The  complexes  were  mostly  supershifted  with 
anti-Jun  or  anti-Fos  antibodies.  Anti-ER  antibody  just  marginally 
decreased  such  complexes,  suggesting  that  a  low  amount  of  ER 
may  be  bound  to  AP-1  sites  under  these  conditions/cells.  More¬ 
over,  DIBA  did  not  display  an  inhibitory  effect  on  E2-  or  TAM- 
stimulated  AP-1 -binding  activity,  possibly  because  ER  binding 


to  DNA  is  not  required  for  its  activity  through  the  nonclassical 
AP-1  pathway  (Jakacka  et  al.,  2001 ;  Webb  et  al.,  1999).  These 
data  further  support  the  specificity  of  DIBA  influencing  ER  bind¬ 
ing  to  DNA. 

DIBA  blocks  occupancy  of  estrogen  target  gene 
promoters  by  ERa 

We  further  used  chromatin  immunoprecipitations  (ChIP)  to  di¬ 
rectly  assess  whether  DIBA  impacts  ERa  binding  to  promoters 
of  estrogen  target  genes.  The  presence  of  the  specific  pro¬ 
moters  in  the  chromatin  immunoprecipitates  was  analyzed  by 
semiquantitative  PCR  by  using  specific  pairs  of  primers  span¬ 
ning  the  estrogen-responsive  regions  in  the  pS2,  c-Myc,  and 
cathepsin  D  (CATD)  gene  promoters  (Figure  3E).  Stimulation 
with  E2  and  TAM  dramatically  increased  ERa’s  occupancy  of 
the  above  three  promoters.  DIBA  remarkably  decreased  such 
occupancy  of  ERa  to  the  target  gene  DNA  sequences  in  chro¬ 
matin.  By  contrast,  ERa  did  not  show  any  interaction  with  the 
distal  promoter  region  (-3351  to  -3551)  of  pS2  promoter. 
These  results  suggested  DIBA  directly  influences  the  ability  of 
ERa  to  bind  to  ERE  in  the  promoter  of  target  genes. 

We  also  used  ChIPs  to  examine  whether  DIBA  affects  ER 
binding  to  AP-1  site  in  a  nontypical  manner  (Figure  3F).  Stim¬ 
ulation  with  E2  and  TAM  induced  a  dramatic  increase  in  the 
occupancy  by  c -fos  or  ERa  of  the  AP-1  site,  but  not  in  the 
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Figure  3.  ER  is  necessary  for  DIBA  to  sensitize  TAM  inhibition 

A:  siRNA-mediated  knockdown  of  ERa  alters  DIBA-mediated  TAM  inhibition  of  resistant  cell  growth.  BT474  cells  were  transduced  with  ERa-siRNA  or  control 
vector  and  incubated  for  96  hr.  Levels  of  ERa  expression  were  examined  by  western  blotting  (inset).  Proliferation  of  the  above  transfected  cells  treated 
with  DIBA  and  4-OH-TAM  in  the  presence  of  E2  was  assayed  by  [3H]thymidine  incorporation.  Data  shown  represent  mean  ±  SEM. 

B-D:  DIBA  inhibits  E2-induced  ERE  (B),  but  not  ARE  (C)  or  AP-1  (D),  DNA  binding.  MCF-7/LCC2  cells  were  treated  with  or  without  5  |iM  DIBA  for  2  hr,  then  stim¬ 
ulated  with  medium  (-),  50  nM  E2,  or  50  nM  4-OH-TAM  (+)  for  20  min.  Nuclear  extracts  were  incubated  in  the  absence  of  antibody,  aER,  aAR,  aJun,  aFos,  or 
normal  rabbit  serum  (NRS)  in  combination  with  32P-labeled  oligonucleotide  probes.  Arrows  indicate  migrational  location  of  each  nonsupershifted  ER,  AR  or 
AP-1  DNA  complex. 

E:  The  recruitment  of  ERa  to  the  promoters  of  estrogen-responsive  genes.  MCF-7/LCC2  cells  were  treated  with  or  without  5  |iM  DIBA  for  2  hr,  then  stimulated  with 
E2, 4-OH-TAM,  or  IGF- 1  for  40  min.  Soluble  chromatin  was  prepared  and  immunoprecipitated  with  anti-ERa.  The  final  DNA  extractions  were  amplified  using  pairs 
of  primers  that  cover  the  regions  of  pS2,  CATD,  and  c-Myc  gene  promoters,  as  indicated.  The  distal  region  (approximately  -3351  to  -3551 )  of  the  pS2  gene 
promoter  was  examined  for  the  presence  of  ERa  (bottom  row). 

F:  The  recruitment  of  ERa  to  the  promoter  of  an  estrogen-induced  AP-1  -dependent  gene  MMP-1 .  Soluble  chromatin  was  immunoprecipitated  with  antibodies 
against  ERa  or  c-Fos.  The  final  DNA  extractions  were  amplified  using  pairs  of  primers  that  cover  the  AP-1  site  as  indicated  or  the  non-AP-1  -specific  site  (approx¬ 
imately  +2555)  of  the  MMP-1  gene  promoter. 
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Figure  4.  DIBA  inhibits  ERE  transactivation 

A-C:  MCF-7/LCC2  orMCF-7  cells  were  transfected  with  a  pGL3-TK-ERE  luciferase,  pGL3-TK-ARE  luciferase,  or  pGL3-AP-l  luciferase  construct,  respectively.  After 
addition  of  4-OH-TAM  (50  nM)  and/or  DIBA  (5  pM)  for  2  hr,  cells  were  stimulated  with  or  without  50  nM  E2  or  100  nM  R1881  for  1 6  hr.  Luciferase  activity  of  lysed 
cells  was  measured  and  normalized.  Data  shown  represent  mean  ±  SEM. 

D:  MDA-MB-468  cells  were  transfected  with  a  wild-type  ER  (pSG5-HE0),  a  series  of  human  ER  deletion  mutants  including  pSG5-HEl  1,  pSG5-HE19,  pSG5-HE16, 
pSG5-HEZF,  or  pSG5  control  plasmids  and  a  pGL3-TK-ERE  luciferase  reporter.  After  24  hr,  the  transfected  cells  were  treated  with  DIBA,  E2, 4-OH-TAM,  and  IGF-1 
for  an  additional  24  hr.  Luciferase  activity  of  lysed  cells  was  measured  and  normalized.  Data  shown  represent  mean  ±  SEM. 


non-AP-1  -specific  site  in  promoter  of  matrix  metalloproteinase 
1  (MM P-1),  an  estrogen-induced/AP-1  -dependent  gene  pro¬ 
moter  containing  AP-1  sites  but  no  ERE  sequences  (DeNardo 
et  al.,  2005).  DIBA  did  not  affect  such  occupancy  of  c -fos  or 
ERa,  consistent  with  the  observation  by  EMSA. 

DIBA  inactivates  ligand-dependent  ERE  transactivation 

To  determine  whether  DIBA  might  affect  TAM-mediated  ER 
transcription  in  TAM-resistant  breast  cancer  cells,  we  tested 
transactivation  of  MCF-7/LCC2  (Figure  4A)  and  MCF-7  (Fig¬ 
ure  4B)  cells  transfected  with  the  ERE-luciferase  reporter 
gene.  E2  activated  ERE  transactivation  in  both  cell  lines.  TAM 
alone  suppressed  E2-induced  ERE  transactivation  in  MCF-7 
cells,  whereas  it  increased  ERE  transactivation  in  MCF-7/ 
LCC2  cells.  DIBA  significantly  reduced  ERE  transactivation 
stimulated  by  4-OH-TAM  and  E2  in  MCF-7/LCC2  cells.  By  con¬ 
trast,  DIBA  did  not  affect  androgen-responsive  element  (ARE) 


transactivation  mediated  by  R1881  in  both  MCF-7  and  MCF- 
7/LCC2  cells  (Figures  4A  and  4B).  Furthermore,  DIBA  did  not 
inhibit  transactivation  of  AP-1-luc  (Figure  4C).  These  data  indi¬ 
cate  that  DIBA  selectively  suppresses  TAM-stimulated  ER 
DNA  binding  and  subsequent  ERE  transactivation. 

To  further  validate  the  target  specificity  of  DIBA  on  ligand- 
dependent  ERE  transcription  in  TAM-resistant  breast  cancer 
cells,  we  cotransfected  the  wild-type  human  ERa  (HE0),  a  series 
of  human  ER  deletion  mutants  (Kumar  et  al.,  1987)  including 
HEZF  (ER  depleted  of  zinc  finger  domains,  ER-AZF),  HE1 1  (ER 
depleted  of  DBD,  ER-ADBD),  HE1 9  (ER  depleted  of  A/B  regions 
but  containing  DBD  and  AF-2  domain,  ER-AA/B),  HE16  (ER 
depleted  of  D/E/F  regions,  ER-AD/E/F),  or  pSG5  control  expres¬ 
sion  plasmid  and  the  ERE-luciferase  reporter  gene  into  the 
ER-negative  MDA-MB-468  cells.  As  shown  in  Figure  4D,  over¬ 
expression  of  ERa,  compared  to  pSG5  control,  remarkably 
resulted  in  ERE  transactivation.  E2  strongly  activated  ERE 
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transactivation,  whereas  TAM  did  not  show  significant  inhibition 
on  such  transactivation.  Also,  induction  of  ERE  transactivation 
by  E2  was  observed  in  the  cells  overexpressed  by  ERa  mutant 
HE19  containing  completed  AF-2  domain  and  DBD.  Deletion 
of  zinc  finger  domains  (HEZF)  or  the  entire  DBD  (HE1 1)  resulted 
in  decreasing  ERE  transactivation  stimulated  by  E2,  suggesting 
that  zinc  fingers  in  DBD  are  required  for  ligand-dependent  ERE 
transactivation.  DIBA  significantly  enhanced  TAM  inhibition  in 
the  wild-type  ERa-  or  a  mutant  HE1 9  (ER-AA/B)-overexpressing 
cells,  but  not  in  HE11  (ER-ADBD)-  or  HEZF  (ER-AZF)-overex- 
pressing  or  control  cells.  Since  DIBA,  as  a  zinc  finger  inhibitor, 
has  been  demonstrated  to  preferentially  disrupt  the  ER  DNA- 
binding  domain,  the  inhibitory  effect  of  DIBA  on  ligand-induced 
ERE  transcription  in  TAM-resistant  breast  cancer  cells  may  be 
related  to  interruption  of  zinc  finger  domains  within  ER-DBD. 

DIBA  reduces  ligand-independent  ERE  transactivation 

The  ligand-independent  ERE  transcription  was  also  measured  in 
the  above  wild-type  ER  or  mutant-transfected  MDA-MB-468 
cells.  As  shown  in  Figure  4D,  in  the  case  of  the  wild-type  ER 
(HEO)-overexpressing  cells,  IGF-1  strongly  induced  ERE  trans¬ 
activation.  Deleting  zinc  finger  domains  or  the  entire  DBD 
decreased  ERE  transactivation  stimulated  by  either  E2  or 
IGF-1,  even  though  this  mutant  contains  completed  AF-1  and 
AF-2  domains.  However,  induction  of  ERE  transactivation  by 
IGF-1  was  observed  in  the  cells  overexpressed  by  the  ER  mu¬ 
tant  HE16  containing  a  completed  N-terminal  A/B  domain  and 
DBD,  suggesting  that  DBD  is  required  for  both  ligand-depen¬ 
dent  and  -independent  ERE  transactivation.  IGF-1 ’s  activation 
of  ERE  transcription  is  not  only  dependent  on  the  AF-1  domain 
itself,  but  is  also  mediated  through  the  interaction  between  DBD 
and  AF-1  domains,  consistent  with  previous  observations  that 
long-range  allosteric  communication  occurs  in  two  separated 
domains  of  the  androgen  receptor  (Brodie  and  McEwan, 
2005),  glucocorticoid  receptor  (Kumar  et  al.,  1999),  and  proges¬ 
terone  receptor  (Bain  et  al.,  2000). 

Moreover,  DIBA  blocked  ERE  transactivation  stimulated  by 
IGF-1  in  the  cells  overexpressed  by  a  wild-type  ER  or  the  ER  mu¬ 
tant  (HE16)  containing  A/B/C  domains.  Such  an  inhibitory  effect 
of  DIBA  was  not  observed  in  the  ER  mutants  HE1 1  (ER-ADBD)- 
and  HEZF  (ER-AZF)-transfected  cells,  indicating  that  inhibition 
of  DIBA  on  the  “steroid-independent  activation”  of  ER  by 
growth  factor  signals  was  related  to  DBD-mediated  intramolec¬ 
ular  communication  with  the  AF-1  domain,  which  may  also  be  in¬ 
volved  in  DIBA  functionally  suppressing  TAM-resistant  breast 
cancer  cells. 

DIBA  decreases  the  TAM-bound  ER  association 
with  AIB1 

Activated  AIB1  probably  translocates  to  nucleus  (Schiff  et  al., 
2004),  where  it  can  interact  with  ER;  therefore,  we  utilized 
a  coimmunoprecipitation  experiment  to  analyze  whether  DIBA 
impacts  the  ERa  interaction  with  AIB1 .  Cell  extracts  were  immu- 
noprecipitated  with  an  anti-ERa-specific  antibody;  immunopre- 
cipitates  were  developed  on  western  blots  with  anti-AIBI  (upper 
panel)  or  anti-ERa  (lower  panel).  In  MCF-7/LCC2  cells  (Fig¬ 
ure  5A),  the  AIB1  can  be  coprecipitated  with  ERa  in  cells  treated 
with  E2,  4-OH-TAM,  or  IGF-1,  indicating  that  a  direct  protein- 
protein  interaction  occurs  between  nuclear  receptor  ERa 
and  its  coactivator  AIB1  upon  addition  of  E2,  4-OH-TAM,  and 
IGF-1.  Notably,  DIBA  significantly  decreased  such  ER 


interaction  with  AIB1.  In  contrast,  E2,  but  not  TAM,  induced 
this  association  between  ERa  and  AIB1  in  MCF-7  cells.  These 
data  support  that  the  effect  of  DIBA  on  TAM-resistant  MCF-7/ 
LCC2  cells  may  be  through  dissociation  of  the  coactivator 
AIB1  complexes  from  TAM-bound  ERa. 

DIBA  increases  association  of  TAM-bound 
ERa  with  NCoR 

Several  lines  of  evidence  indicate  that  the  nuclear  receptor 
corepressor  (NCoR)  complex  mediates  the  inhibitory  effects  of 
TAM  (Keeton  and  Brown,  2005).  Thus,  we  examined  whether 
DIBA  affects  NCoR  modulation  of  the  response  of  ERa  to  TAM 
by  using  a  coimmunoprecipitation  experiment.  Cell  extracts 
were  immunoprecipitated  with  an  anti-NCoR  specific  antibody; 
immunoprecipitates  were  developed  on  western  blots  with  anti- 
ERa  (Figure  5B).  In  control  MCF-7  cells,  TAM  induced  the  asso¬ 
ciation  between  ERa  and  NCoR  while  E2  did  not  affect  it,  which 
may  mediate  the  antagonistic  effect  of  TAM  on  its  sensitive  cells. 
In  MCF-7/LCC2  cells,  a  little  ERa  can  be  coprecipitated  with 
NCoR,  suggesting  that  a  weak  constitutive  interaction  occurs 
between  nuclear  receptor  ERa  and  NCoR,  which  is  consistent 
with  the  previous  observations  that  interactions  of  ERa  with 
NCoR  in  vitro  appear  to  occur  regardless  of  the  ligand  state  of 
the  receptor  (Smith  et  al.,  1997;  Voss  et  al.,  2005).  Although 
E2  and  IGF-1  significantly  decreased  such  interaction,  TAM 
alone  did  not  increase  it.  DIBA  remarkably  increased  NCoR 
association  with  ERa  in  the  presence  of  E2,  IGF-1,  and  TAM, 
suggesting  that  effect  of  DIBA  on  TAM-resistant  MCF-7/LCC2 
cells  may  also  occur  through  association  of  the  corepressor 
NCoR  complexes  with  ERa. 

DIBA  mediates  chromatin-associated  recruitment 
of  ERa  and  cofactors 

To  examine  whether  interaction  between  ERa  and  AIB1  or  ERa 
and  NCoR  is  chromatin  associated,  we  performed  ChIP  assays 
of  ER  followed  by  the  Re-ChIP  analysis  of  either  AIB1  or  NCoR, 
analyzing  the  assembly  of  ERa-cofactor  complex  components 
on  a  well-characterized  estrogen-responsive  pS2  promoter 
(Figure  5C).  The  soluble  chromatin  derived  from  MCF-7/LCC2 
cells  was  subjected  to  ChIP  with  ERa  antibodies;  subsequently, 
the  released  immune  complexes  were  divided  into  two  aliquots 
for  the  Re-ChIP  using  AIB1  antibodies  or  NCoR  antibodies.  The 
same  Re-ChIP  was  also  performed  on  the  unbound  supernatant 
fractions  from  the  primary  immunoprecipitation.  The  ChIP  assay 
of  ERa  antibodies  showed  that  strong  binding  of  ERa  to  the 
pS2  promoter  was  induced  by  E2  or  TAM.  DIBA  significantly 
decreased  E2-  or  TAM-occupied  ERa  binding  to  the  estrogen- 
responsive  DNA  sequences  in  the  pS2  promoter.  The  Re-ChIP 
assay  using  AIB1  antibodies  illustrated  that  E2  or  TAM  induced 
occupancy  of  the  pS2  promoter  by  ER  and  the  coactivator  AIB1 . 
However,  the  Re-ChIP  assay  using  NCoR  antibodies  showed 
that  a  marginal  recruitment  of  the  NCoR  occurred  in  the  absence 
of  ligand,  while  stimulating  E2  or  TAM  abolished  such  promoter 
occupancy  by  ERa-NCoR  complexes,  indicating  that  interac¬ 
tions  between  ERa-AIBI  and  between  ERa-NCoR  are  chroma¬ 
tin  associated.  After  DIBA  treatment,  there  were  very  low  levels 
of  E2-  orTAM-induced  recruitment  of  ERa-AIBI  and  ERa-NcoR 
complexes  to  chromatin.  Combined  with  the  data  obtained  from 
coimmunoprecipitation  experiments  (Figures  5A  and  5B),  these 
results  suggested  that  DIBA-induced  changes  in  ERa  associa¬ 
tion  with  cofactors  led  to  inhibition  of  ERa  binding  to  DNA,  in 
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Figure  5.  Effect  of  DIBA  on  ERa  association  with 
AIB1  or  NCoR 

A  and  B:  ERa  association  with  cofactors  assayed 
by  coimmunoprecipitation.  MCF-7/LCC2  or 
MCF-7  cells  were  treated  with  or  without  DIBA 
for  2  hr,  and  then  stimulated  with  E2,  4-OH-TAM, 
or  IGF-1  for  24  hr  before  lyses.  A:  Western  blotting 
analysis  with  anti-AIBl  (upper  panel)  or  anti-ERa 
(lower  panel)  was  performed  on  anti-ERa  immu- 
noprecipitates.  B:  Western  blotting  analysis  was 
performed  with  anti-ERa  (upper  panel)  or  anti- 
NCoR  (lower  panel)  on  anti-NCoR  immunopre- 
cipitates. 

C:  Recruitment  of  ERa  and  cofactors  assayed  by 
ChIP-Re-ChIP.  Soluble  chromatin  was  immuno- 
precipitated  with  antibodies  against  ERa  (1°  IP). 
The  supernatant  was  collected  and  reimmuno- 
precipitated  with  antibodies  against  AIB1  or 
NCoR  (Supernatant  Re-IP).  Similar  reciprocal 
Re-IPs  were  also  performed  on  complexes  eluted 
from  the  1°  IPs  (Bound  Re-IP). 

D:  Time  course  of  GFP-ERa  redistribution.  MCF-7/ 
LCC2  cells  were  transiently  transfected  with 
pEGFP-C2-hERa.  Live  cells  expressing  GFP-ERa 
were  pretreated  with  vehicle  or  DIBA  for  2  hr, 
followed  by  stimulation  with  50  nM  E2  or  50  nM 
TAM.  Time  courses  of  GFP-ER  distribution  were 
analyzed  at  10  min  intervals.  Scale  bar,  5  |im. 


turn  blocking  transcription  of  target  genes,  which  aided  the 
synergism  between  DIBA  and  TAM. 

Since  cofactor  association  can  influence  ERa  cellular  localiza¬ 
tion,  we  used  a  transcriptionally  active  green  fluorescent  pro- 
tein-ERa  chimera  (GFP-ERa)  to  examine  whether  DIBA  affects 
ERa  cellular  distribution.  MCF-7/LCC2  cells  were  transiently 
transfected  with  pEGFP-C2-hERa,  and  live  cells  expressing 
GFP-ERa  were  analyzed  at  1 0  min  intervals  under  confocal  laser 
scanning  microscopy.  Without  ligand,  GFP-ERa  was  observed 
only  in  the  nucleus,  excluding  the  nucleolus,  with  a  diffuse  dis¬ 
tribution.  Upon  adding  E2,  GFP-ERa  was  dramatically  redistrib¬ 
uted  from  a  reticular  to  punctate  pattern  within  the  nucleus  (Fig¬ 
ure  5D).  A  similar  reorganization  occurred  with  TAM.  In  the  cells 
pretreated  with  DIBA,  neither  E2  nor  TAM  produced  the  above 
apparent  subnuclear  redistribution  patterns.  These  results  dem¬ 
onstrated  that  DIBA  inhibited  E2-  or  TAM-induced  ERa  nuclear 
distribution. 

DIBA  dephosphorylates  ERa  at  serine-167, 
but  not  serine-118 

The  human  ERa  AF-1  function  is  potentiated  by  the  phosphory¬ 
lation  of  serine  residues  of  human  ERa  A/B  domain  after  stimu¬ 
lation  with  its  ligands  and  nonsteroidal  growth  factors  (EGF  and 
IGF-1)  (Lannigan,  2003;  Yamashita  et  al.,  2005).  We  thus  inves¬ 
tigated  whether  DIBA  may  modulate  ERa  phosphorylation  by 


using  site-specific  antiphosphoserine  antibodies  against  ERa 
at  Ser-118  or  Ser-167  (Figure  6A).  E2,  4-OH-TAM,  and  IGF-1 
stimulated  Ser-167  phosphorylation,  whereas  there  was  no  sig¬ 
nificant  difference  in  the  level  of  phosphorylation  of  ERa  at  Ser- 
1 1 8  in  MCF-7/LCC2  cells  with  the  above  treatments.  While  DIBA 
inhibited  phosphorylation  of  ERa  at  Ser-1 67  induced  by  all  stim¬ 
uli,  it  affected  neither  Ser-118  phosphorylation  nor  the  expres¬ 
sion  of  ERa.  It  has  been  demonstrated  that  ERa  phosphorylation 
at  Ser-1 67,  but  not  at  Ser-1 1 8,  conferred  DNA  binding  and  tran¬ 
scriptional  activation  (Joel  et  al.,  1998)  as  wellasTAM  resistance 
(Campbell  et  al.,  2001).  Since  the  structure  of  the  N-terminal 
AF-1  domain  appears  to  be  influenced  by  the  DBD  (Graham 
et  al.,  2000),  and  DIBA  selectively  reacts  with  zinc  finger  of 
ER-DBD,  it  seems  likely  that  DIBA  may  interfere  with  phos¬ 
phorylation  of  Ser-167  in  AF-1  proximal  to  the  DBD  site  through 
intramolecular  communication,  concurring  with  the  above 
observation  on  the  effect  of  DIBA  in  ligand-independent  ERE 
transactivation  (Figure  4D),  and  may  also  contribute  to  DIBA 
sensitizing  the  resistant  cells  to  TAM. 

DIBA  does  not  affect  expression  and  phosphorylation 
of  AIB1  and  MAPK 

Since  ER  coactivator  AIB1 ,  like  ER  itself,  is  phosphorylated  and 
activated  by  different  signaling  kinases,  including  the  p42/44 
MAPK,  which  can  be  activated  by  HER2  (Font  de  Mora  and 
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Figure  6.  Effect  of  DIBA  on  expression  and  phosphorylation  of  ERa,  AIB1 ,  and 
MAPK 

A:  MCF-7/LCC2  cells  were  treated  with  or  without  DIBA  for  2  hr,  then  simu¬ 
lated  with  E2,  4-OH-TAM,  or  IGF-1  for  20  min  before  lyses.  Western  blotting 
analysis  was  performed  with  anti-phospho-ER  (Ser-1 67  or  Ser-1 1 8)  or  anti-ER. 
B:  MCF-7/HER2-18  cell  lysates  were  immunoprecipitated  with  anti-AIBl.  Im- 
munoprecipitates  were  blotted  with  anti-phosphoserine  (upper  panel), 
anti-phosphotyrosine  (middle  panel),  or  anti-AIBl  (lower  panel). 

C  and  D:  Cell  lysates  of  MCF-7/HER2-1 8  (C),  MCF-7/LCC2,  and  BT474  (D)  were 
analyzed  with  anti-phospho-MAPK  for  blot  (upper  panel)  or  anti-MAPK 
(lower  panel)  for  re-blot. 


Brown,  2000),  we  examined  whether  DIBA  attenuates  phos¬ 
phorylation  of  AIB1  in  TAM-resistant  cells  (Figure  6B).  MCF-7/ 
HER2-18  cell  extracts  were  immunoprecipitated  with  an  anti- 
AIBl  specific  antibody;  the  immunoprecipitates  were  developed 
on  western  blots  with  anti-phosphoserine  (upper  panel),  anti- 
phosphotyrosine  (middle  panel),  or  anti-AIBl  (lower  panel). 
The  phosphorylation  of  serine,  but  not  tyrosine,  of  AIB1  could 
be  observed  in  cells  stimulated  with  E2,  4-OH-TAM,  and  EGF. 
However,  DIBA  did  not  affect  such  phosphorylation,  indicating 
that  DIBA  inactivates  ER  Ser-1 67  phosphorylation,  but  does 
not  affect  expression  and  phosphorylation  of  AIB1,  possibly 
due  to  AIB1  ’s  lacking  a  zinc  finger,  although  the  signaling  from 
the  EGFR/HER2  family  activates  ER  and  AIB1  by  the  p42/44 
MAPK. 

We  further  examined  whether  DIBA  disrupted  phosphoryla¬ 
tion  of  MAPK  in  different  TAM-resistant  breast  cancer  cell  lines. 
Figure  6C  shows  the  same  pattern  for  the  phosphorylation  of 
MAPK  as  that  for  AIB1  in  MCF-7/HER2-18  cells.  Similar  results 


Figure  7.  Effect  of  DIBA  on  expression  and  phosphorylation  of  HER2  and  AKT 
A:  MCF-7,  MCF-7/LCC2,  and  MCF-7/HER2-1 8  cells  were  treated  with  DIBA  for 
2  hr  and  then  stimulated  with  E2, 4-OH-TAM,  or  IGF-1  for  20  min.  HER2  expres¬ 
sion  was  analyzed  by  western  blotting  with  anti-HER2. 

B:  MCF-7/LCC2,  MCF-7/HER2-1 8,  and  BT474  cells  were  treated  as  described 
in  A,  except  that  the  antibodies  were  anti-phospho-HER2  for  blot  (upper 
panel)  or  anti-HER2  for  re-blot  (lower  panel). 

C:  MCF-7/LCC2  cells  were  treated  as  described  in  A,  except  that  the  anti¬ 
bodies  were  anti-phospho-AKT  for  blot  or  anti-AKT  for  re-blot. 

were  observed  in  MCF-7/LCC2  and  BT474  cells  (Figure  6D);  the 
ratio  of  phospho-MAPK  to  total  MAPK  was  not  significantly 
changed  after  the  treatment  with  DIBA. 

DIBA  does  not  influence  expression  and  phosphorylation 
of  HER2  and  AKT 

Since  overexpression  of  HER2  and  high  levels  of  phosphory- 
lated  AKT  may  also  contribute  to  TAM  resistance  (Gutierrez 
et  al.,  2005;  Osborne  et  al.,  2003),  we  examined  whether  DIBA 
disrupted  expression  and  phosphorylation  of  HER2  and  AKT 
in  TAM-resistant  breast  cancer  cell  lines.  Compared  to  MCF-7 
cells,  MCF7/HER2-18,  MCF-7/LCC2,  and  BT474  cells  ex¬ 
pressed  a  considerable  level  of  HER2  (Figures  7  A  and  7B).  There 
are  no  significant  changes  in  HER2  expression  in  DIBA-treated 
cells.  Moreover,  even  though  TAM  and  IGF-1  induced  remark¬ 
able  phosphorylation  of  HER2  (Stoica  et  al.,  2000),  DIBA  did 
not  affect  it  (Figure  7B),  nor  did  DIBA  significantly  affect  TAM- 
or  IGF-1 -stimulated  phosphorylation  of  AKT  (Figure  7C)  in 
MCF-7/LCC2  cells.  These  results  indicate  the  inhibitory  effect 
of  DIBA  on  TAM-resistant  cell  proliferation  is  not  based  on  inac¬ 
tivation  of  HER2,  MAPK,  and  PI3-K/AKT. 
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Discussion 

ER  plays  a  major  role  in  many  cases  of  breast  cancer  and  appar¬ 
ently  contributes  to  the  growth  of  TAM-responsive,  acquired 
TAM-resistant,  and  de  novo  ER-positive  resistant  models  (Gee 
et  al.,  2005).  By  using  siRNA  to  deplete  ER  of  BT474,  an  ER-pos¬ 
itive  but  TAM-resistant  breast  cancer  cell  line,  we  found  that 
DIBA  rendered  TAM  inhibition  on  parent  ER-positive  cells,  but 
not  on  ER-depleted  breast  cancer  cells  (Figure  3A),  suggesting 
that  DIBA  function  on  TAM-resistant  breast  cancer  cells  is  ER 
dependent.  We  previously  discovered  that  DIBA  preferentially 
disrupted  the  vulnerable  zinc  fingers  of  the  ERa  DNA-binding 
domain,  thus  blocking  ER  DNA  binding  and  transactivation 
(Wang  et  al.,  2004).  In  this  study,  we  used  chromatin  immuno- 
precipitations  to  directly  assess  ER  binding  to  DNA  on  estrogen 
target  genes.  Stimulation  with  E2  and  TAM  dramatically  in¬ 
creased  ER’s  occupancy  of  the  pS2,  c-Myc,  and  CATD  gene 
promoters.  DIBA  remarkably  decreased  such  occupancy  of 
ER  to  the  target  DNA  sequences  in  chromatin  (Figure  3E).  These 
results  suggested  DIBA  directly  influences  the  ability  of  ER  to 
bind  to  DNA,  consistent  with  the  data  obtained  from  EMSA 
(Figure  3B).  Moreover,  DIBA  resulted  in  inhibition  of  ligand- 
dependent  and  -independent  ERE  transactivation  (Figure  4D). 
Therefore,  targeted  disruption  of  ER  is  necessary  for  DIBA, 
a  zinc  finger  inhibitor,  to  sensitize  TAM  inhibition  of  resistant 
breast  cancer  cells  through  interfering  with  ER  DNA  binding 
and  subsequent  ERE  transactivation. 

Nuclear  receptor  function  is  modulated  by  transcriptional  cor¬ 
egulators  (Klinge  et  al.,  2004;  McKenna  and  O’Malley,  2002; 
Shang  et  al.,  2000;  Shang  and  Brown,  2002;  Tikkanen  et  al., 
2000).  The  relative  level  of  these  coactivators  and  corepressors 
might  determine  the  balance  of  agonist  and  antagonist  proper¬ 
ties  of  TAM.  Here,  we  used  coimmunoprecipitation  to  clarify  that 
DIBA  decreased  physical  association  of  TAM-bound  ER  with  its 
coactivator  AIB1  (Figure  5A),  whereas  it  increased  ER  interac¬ 
tion  with  its  corepressor  NCoR  (Figure  5B).  Moreover,  ChIP  ex¬ 
periments  of  ERa  followed  by  either  Re-ChIP  of  AIB1  or  NCoR 
also  showed  that  effect  of  DIBA  on  E2-  or  TAM-induced  associ¬ 
ation  between  ER  and  AIB1  and  dissociation  between  ERa  and 
NCoR  in  the  level  of  chromatin  (Figure  5C),  suggesting  that 
Dl BA-mediated  changes  in  ERa  interaction  with  cofactors 
resulted  in  blockage  of  TAM-bound  ERa  binding  to  targeted 
gene  promoter  and  transcription.  Notably,  the  ER-cofactor 
association  caused  by  DIBA  further  influences  TAM-bound  ER 
nuclear  distribution  (Figure  5D),  indicating  other  functional 
changes  of  ERa  may  have  with  chromatin  on/off  rates  or  shut¬ 
tling.  The  above  molecular  mechanisms,  by  which  the  syner¬ 
gism  between  DIBA  and  TAM  impacted  ER  activity,  contributed 
to  DIBA  restoring  TAM’s  antagonist  action  on  TAM-resistant 
breast  cancer  cells  (Figures  1  and  2).  It  may  be  important  to 
note  that  in  our  previous  report,  the  ED50  for  DIBA  ejection  of 
zinc  from  recombinant  ERa  was  25  jiM  (Wang  et  al.,  2004). 
Here  we  show  effects  on  disrupting  ERa/AIBI  or  enhancing 
ERa/NCoR  functions  at  5-fold  less,  5  jiM,  suggesting  a  range 
of  molecular  effects  on  ERa  functions. 

Several  peptide  growth  factors  and  their  intracellular  signaling 
kinases,  notably  MAPK  and  AKT  (Albanell  and  Baselga,  2001), 
have  been  shown  to  mediate  cell  proliferative  responses  and 
phosphorylate  ERa  on  various  AF-1  residues,  promoting  ERa 
transcriptional  activity  in  a  ligand-independent  manner  (Martin 
et  al.,  2000).  In  the  case  of  TAM-resistant  cells,  we  observed 


that  exogenous  IGF-1  stimulated  phosphorylation  of  MAPK  and 
AKT  as  well  as  ERa.  Although  DIBA  did  not  affect  phos¬ 
phorylation  of  ERa  at  Ser-118,  MAPK,  or  AKT,  DIBA  markedly 
inhibited  phosphorylation  of  ERa  at  Ser-167  (Figure  6A),  sug¬ 
gesting  that  inhibitory  effects  of  DIBA  on  a  powerful  functional 
crosstalk  engaged  by  the  IGF-1  and  ER  pathways  may  occur 
through  dephosphorylating  ERa  at  Ser-167.  Thus,  DIBA  disrup¬ 
tion  of  ER  zinc  fingers  resulted  in  not  only  perturbing  DBD- 
dependent  ERE  transactivation  (Figures  4A-4D),  but  also  inter¬ 
fering  with  intramolecular  communication  between  DBD  and 
the  N-terminal  AF-1  domain  (Figure  4D)  to  downregulate  phos¬ 
phorylation  of  Ser-167  (Figure  6A)  induced  by  nonestrogenic 
stimulation  in  TAM-resistant  breast  carcinoma  cells. 

Overexpression  of  HER2  and  high  levels  of  phosphorylated 
AKT  or  ERK1/2  MAPK  may  also  contribute  to  TAM  resistance. 
MCF-7  cells  stably  transfected  with  HER2  (MCF-7/HER2-18) 
are  de  novo  resistant  to  TAM,  in  contrast  to  their  low-expressing, 
responsive  MCF-7  counterparts  (Benz  et  al.,  1993;  Konecny 
et  al.,  2003;  Shou  et  al.,  2004).  Importantly,  EGFR/HER2  signal¬ 
ing  remains  dependent  on  the  ER  in  MCF-7/HER2-18  cells,  as 
evidenced  by  their  retained  sensitivity  to  estrogen  deprivation 
(Shou  et  al.,  2004).  DIBA  did  not  inhibit  phosphorylation  of 
HER2,  MAPK,  and  AKT  (Figures  6  and  7)  or  of  the  coactivator 
AIB1  (Figure  6B).  Moreover,  DIBA  did  not  display  suppression 
of  estrogen-  or  TAM-induced  DNA  binding  (Figures  3D  and  3F) 
and  transactivation  for  AP-1  (Figures  3F  and  4C),  a  nontypical 
ER-binding  site.  These  results  suggested  that  nongenomic  ac¬ 
tions  of  ERa  may  be  not  involved  in  the  synergism  between 
DIBA  and  TAM.  However,  this  idea  can  not  be  totally  excluded; 
most  of  the  data  suggest  that  DIBA  blocks  classical  genomic 
sites  of  ERa. 

In  conclusion,  DIBA  resulted  not  only  in  inhibition  of  ligand-de- 
pendent  ERa  DNA  binding  and  transcription,  but  also  in  effects 
on  ligand-independent  ERE  transactivation.  Of  particular  impor¬ 
tance  was  the  synergism  between  DIBA  and  TAM  in  regulating  re¬ 
cruitment  of  cofactors  to  chromatin  (decreasing  the  interaction  of 
ERa  with  AIB1  and  blocking  dissociation  between  ERa  and  NCoR 
caused  by  E2  or  TAM).  Consequently,  DIBA  restores  the  antago¬ 
nistic  action  of  TAM  in  breast  cancer  cells  that  have  acquired  re¬ 
sistance,  in  turn  quenching  target  gene  expression  and  blocking 
cell  growth  of  TAM-resistant  breast  cancer  cells.  These  studies 
suggest  a  possible  new  approach  in  modifying  TAM  resistance 
and  a  potential  role  for  small  electrophilic  compounds  that  can 
modify  the  particularly  vulnerable  zinc  finger  in  ERa. 

Experimental  procedures 
Cell  and  cell  culture 

The  electrophilic  compound  DIBA  (NSC654077)  was  from  the  Laboratory  of 
Cell  Biology,  National  Cancer  Institute.  The  human  breast  carcinoma  cell  lines 
MCF-7,  ZR-75,  and  MDA-MB-468  were  obtained  from  ATCC  (Manassas,  VA). 
The  MCF-7/LCC2  cell  line  was  from  Dr.  R.  Clarke.  MCF-7/HER2-18  and 
BT474  cell  lines  were  from  Dr.  K.  Osborne.  4-OFI-TAM  and  17(3-Estrodial 
were  purchased  from  Sigma-Aldrich  (St  Louis,  MO).  IC1 182780  was  from  Toc- 
ris  (Ellisville,  MO).  In  experiments  with  estrogen  or  TAM,  cells  were  cultured  in 
phenol  red-free  and  DMEM  or  RPMI  1640  supplemented  with  5%  charcoal- 
dextran-stripped  fetal  calf  serum  for  at  least  2  days. 

Cell  proliferation  and  cell-cycle  analysis 

Cell  proliferation  was  examined  by  measuring  DNA  synthesis  with  [3H]thymi- 
dine  uptake  (Wang  et  al.,  2004).  Cell  cycle  was  analyzed  by  propidium  iodide 
staining  and  FACS  (Li  et  al.,  2006). 
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Electrophoretic  mobility  shift  assay 

Electrophoretic  mobility  shift  assay  (EMSA)  was  performed  as  described 
previously  (Wang  et  al.,  2003).  End-labeled  [32P]  oligonucleotide  probes 
correspond  to  the  ERE  consensus  sequence:  S'-GATCCGTCAGGTCAC 
AGTGACCTGATGGATC-3',  ARE  consensus:  5'-GAAGTCTGGTACAGG 
GTGTTCTTT1  I G-3',  and  AP-1  consensus:  S'-CGCTTGATGAGTCAGCCG 
GAA-3',  respectively. 

Expression  plasmids 

pSG5-HE0,  pSG5-HE11,  pSG5-HE16,  or  pSG5-HE19  expression  plasmids 
were  kindly  provided  by  Dr.  P.  Chambon,  Universite  Louis  Pasteur,  France. 
pSG5-HEZF  was  created  by  site-directed  mutagenesis  of  FIE0  using  the  oli- 
gonuecleotide:  3305:  S'-CTCACTATAGGGCGAATTCCGGCCACGGACCAT 
G ACCAT G  ACCC-3' ;  3306:  5' -CAT AT AGT CGTTAT GTCCTT G AAT ACTT CT C 
TT G AAG AAG G CCTT GTAG C G AGT CTCCTT G G CAG ATTCC-3' ;  3307:  5'-GG 
CCTTCTTCAAGAGAAGTATTCAAGGACATAACGACTATATGTACGAAGTGG 
G  AAT  GAT  G  AAAG  GTG  G  G  -3' ;  3308:  5'  -T  CAGACT  GTGGCAGG  GAAACCC 
TCTGCCTCCCCC-3';  3309:  5'  -  AACT  C  G  AG  CTG  G  ATCCT  CAGACT  GTG  G  C 
AGGGAAACCCTCTGCCTCCCCC-3'  resulting  in  deletion  of  amino  acids 
185-205  and  221-245. 

Transfection  of  SiRNA  for  ERa 

BT474  cells  were  transfected  with  an  ERa-SiRNA  construct  or  control  vector 
for  96  hr  according  to  the  manufacturer’s  instructions  (New  England  BioLabs, 
MA).  Efficacy  of  the  constructs  was  tested  through  western  blot  analysis  of 
the  respective  target  ERa  in  transfected  cells. 

Transfection  of  luciferase  reporter  plasmids 

FuGene-6  was  used  for  transfection  of  luciferase  reporter  plasmids  or 
cotransfection  of  reporter  gene  plasmids  with  ER  expression  plasmids.  Lucif¬ 
erase  assays  were  performed  according  to  the  manufacturer’s  instructions 
(BD  PharM ingen,  San  Diego,  CA). 

Coimmunoprecipitation  and  western  blot  analysis 

Coimmunoprecipitation  and  western  blots  were  performed  as  previously 
described  (Yang  et  al.,  2000).  Antibodies  against  ER,  phospho-ER,  AIB1, 
phospho-AIBI ,  NCoR,  HER2,  phospho-HER2,  AKT,  phospho-AKT,  MAPK, 
phospho-MAPK,  phosphotyrosine,  and  phosphoserine  were  from  Upstate 
Biotechnology  (Lake  Placid,  NY). 

Chromatin  immunoprecipitation 

The  ChIP  assays  were  based  on  a  protocol  described  by  Shang  et  al.  (2000). 
Cells  were  fixed  by  formaldehyde.  Purified  chromatin  samples  were  immuno- 
precipitated  with  anti-ERa  antibody.  DNA,  isolated  from  immunoprecipitated 
material  following  reversal  of  formaldehyde  crosslinking,  was  amplified  by 
PCR.  Promoter-specific  primers  included:  pS2,  5'-CCGGCCATCTCTCAC 
TAT-3'  (forward  primer)  and  S'-ATCTTGGCTGAGGGATCT-S'  (reverse 
primer);  pS2  upstream  primer  pair  for  negative  control,  5'-GAAGACTCCG 
CACCT CAG AC-3'  (forward  primer)  and  S'-CCCTTGTGGGGAATCTGG-S' 
(reverse  primer);  c-Myc,  S'-CCGCCTGCGATGATTTATAC-S'  (forward 
primer)  and  5'-AAGGTGGGGAGGAGACTCAG-3'  (reverse  primer);  Cathep- 
sin  D,  5'  -T  CC  AG  AC  ATCCT  CT  CT  G  GAA-3'  (forward  primer),  5'-GGAGCGG 
AGGGTCCATTC-3'  (reverse  primer).  MMP-1  promote,  5'-TTGCAACACCAA 
GTGATTCCA-3'  (forward  primer)  and  S'-CCCAGCCTCTTGCTACTCCA-S' 
(reverse  primer);  MMP-1  non-AP-1  specific  site,  S'-GAGTACAACTTACA 
TCGTGTTGCAG-3'  (forward  primer)  and  5'-ATATGGCTTGGATGCCATCA 
ATGTC-3'  (forward  primer). 

ChIP  Re-ChIP 

Complexes  were  eluted  from  the  primary  immunoprecipitation  by  incubation 
with  10  mM  DTT  at  37°C  for  30  min  and  diluted  1:50  in  buffer  (1%  Triton 
X-100,  2  mM  EDTA,  150  mM  NaCI,  20  mM  Tris-HCI  [pH  8.1]),  followed  by 
reimmunoprecipitation  with  the  second  antibodies  (Shang  et  al.,  2000). 
ChIP  Re-ChIPs  of  supernatants  were  done  essentially  as  the  primary  IPs. 

Live  microscopy 

MCF-7/LCC2  cells  were  grown  on  14  mm  coverslips  in  35  mm  plates  and 
transfected  with  a  pEGFP-C2-hERa  construct  using  FuGene-6.  Before  ligand 
addition,  the  starved  cells  were  pretreated  with  DIBA  for  2  hr.  Images  were 


acquired  at  1 0  min  intervals  with  a  Zeiss  LSM  51 0  confocal  microscope  using 
a  40x/1 .3  NA  oil  immersion  objective  lens  (Stenoien  et  al.,  2000). 

Human  tumor  xenografts 

Human  MCF-7/LCC2-derived  tumor  xenografts  were  established  in  female 
athymic  Ncr-nu/nu  nude  mice  (National  Cancer  Institute,  Frederick,  MD)  as 
described  previously  (Brunner  et  al.,  1993;  Wang  et  al.,  2004).  Tumor  volume 
is  calculated  as  a2  x  b  x  0.5,  where  “a”  is  the  width  and  “b”  is  the  length  of 
the  tumor.  Formalin-fixed  tissue  sections  were  embedded  in  paraffin,  stained 
with  hematoxylin  and  eosin,  and  examined  under  a  light  microscope.  Animal 
experimentation  was  reviewed  and  approved  by  NCI’s  Animal  Research 
Committee. 
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Abstract  Loss  or  mutation  of  known  tumor  suppressor  genes  accounts  for  a  small  proportion  of  all  breast 

cancers.  We  have  recently  shown  that  interferon  regulatory  factor  1  ( IRF1 )  is  a  putative  tumor  sup¬ 
pressor  gene  in  breast  cancer.  We  now  report  that  the  A4396G  single  nucleotide  polymorphism  in 
the  IRF1  gene  is  more  frequent  in  human  breast  cancer  cell  lines  than  in  the  general  population 
( P  =  0.01).  Furthermore,  A4396G  is  more  frequently  expressed  in  African  American  (black)  than 
in  European  ancestry  (white)  subjects  ( n  =  70  subjects;  P  ^  0.001),  leading  to  a  significant  differ¬ 
ence  in  genotype  distribution  between  these  populations  (P  =  0.002).  ©  2007  Elsevier  Inc.  All 

rights  reserved. 


1.  Introduction 

The  precise  molecular  events  responsible  for  affecting 
breast  cancer  risk  and  disease  progression  remain  to  be  es¬ 
tablished.  Although  an  increasing  number  of  oncogenes 
have  been  identified,  relatively  few  tumor  suppressor  genes 
have  been  implicated  in  driving  the  development  and  pro¬ 
gression  of  this  disease.  The  transcription  factor  interferon 
regulatory  factor  1  gene  (HUGO  symbol  IRFl ;  http:// 
www.gene.ucl.ac.uk/nomenclature)  is  lost,  mutated,  or  rear¬ 
ranged  in  several  cancers  including  some  hematopoietic  [1] 
and  gastric  cancers  [2] .  Several  single  nucleotide  polymor¬ 
phisms  have  also  been  reported  [3,4].  We  have  previously 
shown  that  IRFl  is  associated  with  acquired  estrogen- 
independence  [5]  and  with  antiestrogen  resistance  in  breast 
cancer  [6];  others  have  shown  its  importance  in  normal 
mammary  epithelial  cells  [7]  and  in  other  breast  cancer 
models  [8].  The  IRFl  protein  is  readily  detected  in  breast 
tumors  [9,10]  and  is  expressed  in  a  pattern  consistent  with 
a  putative  gene  network  that  regulates  endocrine  respon¬ 
siveness  in  breast  cancer  cells  [10,11]. 

We  have  shown  that  IRFl  acts  as  a  breast  cancer  sup¬ 
pressor  gene,  repressing  both  the  growth  of  human  breast 
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cancer  cell  xenografts  in  athymic  nude  mice  and  cell  prolif¬ 
eration  in  vitro  [12].  These  tumor  suppressor  and  antiproli¬ 
ferative  activities  are  associated  with  the  ability  of  IRFl  to 
modulate  apoptosis  through  regulation  of  a  caspase  cascade 
[7,12].  Loss  of  heterozygosity  at  5q31.1,  which  includes  the 
IRFl  gene  locus,  has  recently  been  reported  in  both  spo¬ 
radic  and  inherited  breast  cancers  [13,14]. 

We  now  report  that  a  polymorphism  in  the  IRFl  gene,  an 
A— >G  single  nucleotide  polymorphism  at  base  pair  4396 
(A4396G),  is  detected  at  a  higher  frequency  in  human 
breast  cancer  cell  lines  than  in  the  normal  population  and 
is  more  frequently  expressed  in  African  Americans  than 
in  European-ancestry  whites. 


2.  Materials  and  methods 

MCF-7  cells  were  originally  obtained  from  Dr.  Marvin 
Rich  (Michigan  Cancer  Foundation,  Detroit,  MI).  T47D, 
ZR-75-1,  MDA-MB-231,  MDA-MB-435,  and  A1N4  cells 
were  obtained  from  the  Tissue  Culture  Shared  Resource 
of  the  Georgetown  University  Lombardi  Comprehensive 
Cancer  Center.  All  other  cell  lines  were  obtained  from 
the  American  Type  Culture  Collection  (ATCC,  Manassas, 
VA).  Unless  otherwise  indicated,  cell  lines  were  routinely 
grown  in  improved  minimal  essential  medium  (IMEM;  Bi¬ 
ofluids,  Rockville,  MD)  with  phenol  red  and  supplemented 
with  5%  fetal  bovine  serum  (Gibco  Life  Technologies— 
Invitrogen,  Carlsbad,  CA).  A1N4  cells  were  grown  in 
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68  IMEM  supplemented  with  0.5%  fetal  calf  serum,  0.5  pg/ 

69  mL  hydrocortisone,  5  pg/mL  insulin,  and  10  ng/mL  epider- 
70QS|  mal  growth  factor  [15].  All  cells  were  maintained  in  a  hu- 

71  midified  incubator  at  37°C  in  an  atmosphere  containing 

72  95%  air,  5%  C02. 

73  Restriction  fragment  length  polymorphism  (RFLP)  anal- 

74  ysis  was  performed  on  a  fragment  of  the  IRF1  gene  ampli- 

75  tied  from  genomic  DNA  by  polymerase  chain  reaction 

76  (PCR).  For  DNA  extraction,  cells  were  plated  in  T-75cm2 

77  plastic  tissue  culture  flasks  at  a  density  of  1  x  106  cells/ 

78  flask  and  grown  for  24  hours  prior  to  DNA  isolation. 

79  DNA  was  extracted  from  proliferating  subconfluent  cells 
[Q2]  using  the  TRIzol  reagent  (Fife  Technologies,  Gaithersburg, 

80  and  quantified  by  comparing  the  optical  density  ratios 

8 1  (OD26o/OD28o)  obtained  spectrophotometrically  using 

82  a  Beckman  DU640  spectrophotometer  (Beckman,  Fuller- 

83  ton,  CA).  For  RFFP  analysis,  IRF1  specific  primers  were 

84  used  to  amplify  a  portion  of  exon  7  in  the  IRF1  gene  from 

85  genomic  DNA.  The  IRF1  primers  were  bp-4358  (sense), 

86  5/-TTGACCTGTGGCTTCTGCTGT-3/,  and  bp-4639  (anti- 

87  sense),  5/-GTCGCTTGCCTCCCCCTATG-3/  from  Gen- 

88  Bank  (accession  no.  F05072;  http://www.ncbi.nlm.nih. 

89  gov);  100  ng  of  the  appropriate  genomic  DNA  was  used 

90  as  template.  PCR  conditions  were  1  cycle  for  3  min  at 

91  94°C,  35  cycles  of  1  min  at  94°C  and  1  min  at  59°C,  and 

92  1  min  45  sec  at  72°C.  Purified  PCR  product  was  digested 

93  with  the  restriction  endonuclease  Nlalll  (New  England 

94  BioFabs,  Beverly,  MA),  size  fractionated  on  a  20%  Tris— 

95  boric  acid— EDTA  polyacrylamide  gel,  and  stained  with 

96  ethidium  bromide.  Nlalll  is  site- specific  (recognition  se- 

97  quence:  CATG)  and  will  not  cleave  variant  DNA  at  the 

98  4396  base  pair.  Thus,  two  bands  correspond  to  a  homo- 

99  zygote  A— >G  polymorphism  at  4396,  three  bands  signify 

100  a  homozygous  wild  type,  and  four  bands  indicate 

101  a  heterozygote. 

102  For  statistical  analyses,  the  proportions  of  genotypes  and 

103  allele  frequencies  in  cell  lines  and  human  subjects  were 

104  compared  by  y2  analysis.  Unless  otherwise  indicated,  all 

105  probabilities  are  two-tailed;  the  conventional  assessment 

106  ofP  <  0.05  was  applied  to  establish  statistical  significance. 

107 

108 
109 

1  3.  Results  and  discussion 

1 1 1  IRF1  is  implicated  as  a  tumor  suppressor  gene  in  breast 

112  cancer  [12]  and  in  several  hematopoietic  [1]  and  gastric 

113  cancers  [2].  We  found  the  A4396G  polymorphism  in  the 

114  IRF1  gene  when  sequencing  PCR  products  from  MCF-7 

115  human  breast  cancer  cells.  Using  RFFP  analysis,  we  mea- 

116  sured  the  prevalence  of  the  polymorphism  in  breast  cancer 

117  cell  lines  (Table  1;  n  =  17),  other  cancer  cell  lines  (Table  2; 

118  n  =  40),  normal  cell  lines  (Tables  1  and  2;  n  =  5),  and  in 

119  DNA  from  normal  volunteers  obtained  from  the  Coriell 

120  DNA  repository  (Coriell  Institute  for  Medical  Resarch, 

121  Camden,  NJ):  34  African  Americans  and  36  whites  ( n  = 

122  70)  (Table  3).  Allele  frequencies  are  significantly  different 


Table  1 

A4396G  polymorphism  in  breast  cancer  and  normal  mammary  epithelial 
cell  lines 


Cell  line 

bp  4396 

Ethnicitya 

MCF7 

G/G 

Eur 

MCF7 

A/G 

Eur 

NIH/ADR-RES 

A/A 

unknown 

ZR-75-1 

A/A 

Eur 

ML-20  (MCF-7  transfected  with  syk) 

A/G 

Eur 

MKL-4  (MCF-7  transfected  with  FGF4) 

A/G 

Eur 

BRC-230 

A/A 

unknown 

BT-474 

A/G 

Eur 

BT-483 

A/A 

Eur 

BT-549 

G/G 

Eur 

DU4475 

A/A 

Eur 

Hs  578T 

A/A 

Eur 

MDA-MB-134-VI 

A/A 

Eur 

MDA-MB-157 

G/G 

AA 

MDA-MB-231 

G/G 

Eur 

MDA-MB-415 

A/A 

Eur 

MDA-MB-330 

G/G 

Eur 

MDA-MB-453 

A/A 

Eur 

MDA435/LCC6 

G/G 

unknown 

Normal  mammary  cell  lines 

A1N4 

A/A 

MCF10A 

A/A 

Abbreviations:  A  A,  African  American  (black);  Eur,  European  ancestry 
(white). 

a  Ethnicity  data  were  obtained  from  the  original  publications  or  as 
described  by  ATCC. 


between  African  Americans  (y2  test,  P  =  0.004;  higher  fre¬ 
quency  of  G  vs.  A)  and  whites  (y2  test,  P  <  0.001;  higher 
frequency  of  A  vs.  G).  The  distribution  of  genotypes  (allelic 
frequency  of  A  was  0.37  in  African  Americans  and  0.65  in 
whites)  was  significantly  different  between  African  Ameri¬ 
cans  and  whites  ( y 2  test,  P  =  0.002). 

We  compared  allele  frequencies  among  the  breast  cancer 
cell  lines  for  which  we  could  identify  ethnic  origin  as  white 
from  either  the  ATCC  (Manassas,  VA)  or  the  original 
publications  (Table  1)  with  the  white  population  data  from 
Table  3.  We  excluded  MDA-MB-157  cells  obtained  from 
an  African  American  woman,  the  MCF-7  cells  because 
we  found  two  different  genotypes,  and  NIH/ADR-RES 
cells  because  the  precise  origin  of  these  cells  is  uncertain 
[16].  We  excluded  MDA435/FCC6  cells  [17];  the  breast 
origin  of  the  parental  MDA-MB-435  cell  line  has  been 
questioned  [18],  and  in  a  recent  study  we  could  not  show 
similarity  in  the  transcrip  tomes  of  MDA-MB-435  cells 
and  a  series  of  breast  cancer  biopsies,  whereas  other  breast 
cancer  cell  lines  show  significant  similarities  with  these 
breast  tumors  [19].  These  cells  now  appear  to  be  defini¬ 
tively  of  melanoma  origin  [20].  We  also  excluded  two  nor¬ 
mal  mammary  epithelial  cell  lines,  A1N4  [21]  and 
MCF10A  [22].  The  remaining  breast  cancer  cell  lines 
(n  =  11)  have  a  genotype  proportion  significantly  different 
from  that  of  the  white  population  ( y 2  test;  P  =  0.01);  geno¬ 
type  frequencies  in  the  breast  cancer  cell  lines  are  not 
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179  Table  2 

i  oq  A4396G  polymorphism  in  cancer  and  other  cell  lines  not  of  mammary 


[Q6]  origin 


182 

Cell  line 

Tissue  of  origin 

bp  4396 

183 

CCF  STTG1 

Astrocytoma 

A/A 

184 

HeLa 

Cervical  adenocarcinoma 

A/G 

1 

LS  147T 

Colorectal  adenocarcinoma 

A/G 

lOJ 

LS  180 

Colon  adenocarcinoma 

A/G 

186 

CaCo-2 

Colon  adenocarcinoma 

A/A 

187 

HCT-116 

Colon  adenocarcinoma 

A/A 

188 

HCT-15 

Colorectal  adenocarcinoma 

G/G 

189 

HEC-l-A 

Endometrial  adenocarcinoma 

A/A 

190 

A-431 

Epidermoid  squamous  carcinoma 

A/A 

Hs  913T 

Fibrosarcoma 

A/G 

191 

KATO  III 

Gastric  adenocarcinoma 

A/G 

192 

A- 172 

Glioblastoma 

A/A 

193 

HepG2 

Hepatocellular  carcinoma 

A/A 

194 

K-562 

Leukemia  (chronic  myeloid  blast  crisis) 

A/A 

195 

Jurkat 

Leukemia  (T-Lymphocyte) 

G/G 

IMR-90 

Lung  fibroblasts  (normal) 

A/G 

196 

A-549 

Lung  (bronchoalveolar  carcinoma) 

A/A 

197 

Calu-3 

Lung  adenocarcinoma 

A/A 

198 

NCI-H209 

Lung  adenocarcinoma 

A/A 

199 

NCI-H345 

Lung  carcinoma 

A/A 

200 

NCI-H520 

Lung  squamous  cells  carcinoma 

G/G 

MOLT-4 

Lymphoblastic  leukemia  (T-cell) 

A/A 

201 

MOLT- 3 

Lymphoblastic  leukemia  (T-cell) 

A/A 

202 

CCRF-CEM 

Lymphoblastic  leukemia  (T-cell) 

A/A 

203 

CCRF-HSB-2 

Lymphoblastic  leukemia  (T-cell) 

A/A 

204 

CCRF-SB 

Lymphoblastic  leukemia  (T-cell) 

A/A 

205 

Daudi 

Lymphoma  (Burkitt) 

G/G 

NHL 

Lymphoma  (Non-Hodgkin) 

A/G 

206 

IMR-32 

Neuroblastoma 

A/G 

207 

BE  (2)  M17 

Neuroblastoma 

G/G 

208 

Hs  683 

Neuroglioma 

A/G 

209 

CaOV-3 

Ovarian  adenocarcinoma 

G/G 

91  0 

COLO-357 

Pancreatic  adenocarcinoma 

G/G 

Z1U 

FaDu 

Pharyngeal  carcinoma 

A/A 

211 

Be  Wo 

Choriocarcinoma 

A/G 

212 

JAR 

Choriocarcinoma 

A/G 

213 

JEG-3 

Choriocarcinoma 

A/G 

214 

HL-60 

Promyelocytic  leukemia 

A/A 

91  S 

LNCaP 

Prostate  adenocarcinoma 

G/G 

Z 1  J 

DU  145 

Prostate  carcinoma 

A/A 

216 

A-204 

Rhabdomyosarcoma 

A/A 

217 

A-673 

Rhabdomyosarcoma 

G/G 

218 

Hs-27 

Foreskin  fibroblast  (normal) 

A/A 

219 

Hs-68 

Foreskin  fibroblast  (normal) 

A/G 

220 

221 

222 

223 

significantly  different  from  those  in  a  random  sample  of 

224 

other  cancer  cell  lines  (Table  2;  n  =  40).  Of  the  excluded 

225 

breast  cancer  cell  lines,  only  the  NIH/ADR-RES  genotype 

226 

is  wild  type,  and  the  outcome  is  unaffected  when  all  the 

227 

cancer  cell  lines  from  Table  2  are  compared. 

228 

229 

230 

Table  3 

^  1  [Q7]  A4396G  polymorphism  in  healthy  human  volunteers 

232 

Ethnicity 

A/A  A/G 

G/G 

233 

African  American  (n  =  34)  3  19 

12 

234 

European  ancestry  (n  =  36)  14  19 

3 

The  greater  prevalence  of  A4396G  in  breast  cancer  cell  235 

lines  derived  from  whites,  compared  with  that  in  the  normal  236 

white  population,  suggests  an  association  with  breast  cancer.  237 

The  genotype  proportions  in  breast  cancer  cells  are  also  seen  238 

in  other  cancer  cell  lines,  further  suggesting  an  association  of  239 

A4396G  with  cancer.  Nonetheless,  we  cannot  fully  exclude  240 

the  possibility  that  this  is  a  cell  culture  artifact.  Our  identifi-  241 

cation  of  two  MCF-7  genotypes  (A/G;  G/G)  suggests  genetic  242 

drift  in  vitro,  which  may  contribute  to  the  phenotypic  diver-  243 

sity  of  this  cell  line.  Ongoing  studies  are  now  examining  244 

the  prevalence  of  this  polymorphism  in  breast  tumors.  245 

How  the  A4396G  polymorphism  contributes  to  the  tumor  246 

suppressor  role  of  IRF1  in  breast  cancer  is  unknown.  Noguchi  247 

et  al.  [4]  first  identified  the  A4396G  polymorphism  in  periph-  248 

eral  blood  lymphocytes  from  several  atopic  and  asthmatic  249 

families  of  Japanese  descent,  although  neither  allele  was  sig-  250 

nificantly  associated  with  transmission  to  asthmatic  children.  25 1 

The  switch  in  nucleotide  use  does  not  change  the  amino  acid  252 

sequence  of  the  translated  protein.  Seven  splice  variants  have  253 

been  previously  reported  in  the  IRF1  gene  [23];  alternative  254 

splicing  in  exons  7,8,  and  9  negatively  regulates  IRF1  in  cer-  255 

vical  cancer,  and  this  splicing  likely  affects  its  tumor  suppres-  256 

sor  activities  [24].  It  is  not  known  if  A4396G  is  an  active  257 

splice  site  but  in  silico  analysis  using  GeneSplicer  [25]  sug-  258 

gests  that  this  polymorphism  is  unlikely  to  generate  a  novel  259 

splice  site.  In  contrast,  the  A4396G  polymorphism  may  af-  260 

feet  putative  transcription  factor  binding  at  internal  sites  261 

[26] .  Compared  with  the  wild  type  allele,  A4396G  loses  bind-  262 

ing  sites  for  microphthalmia  transcription  factor  ( MITFI  263 

TFE3 ),  paraxis  (TCF15),  neurogenin  1  and  3,  and  Myc-  264 

Max  heterodimers;  these  sites  are  replaced  with  a  single  hairy  265 

and  enhancer  of  split  1  ( HES1 )  site.  In  melanoma,  MITF  tran-  266 

scriptionally  activates  the  tumor  suppressor  and  cell  cycle  in-  267 

hibitor  INK4A,  leading  to  cell  cycle  arrest  [27].  Whether  268 

MITF  or  any  of  these  other  transcription  factors  differentially  269 

regulates  the  wild  type  versus  A4396G IRF1  allele  is  under  270 

investigation,  because  altered  regulation  of  IRF1  tumor  sup-  27 1 

pressor  activities  could  significantly  affect  breast  cancer  risk.  272 

Of  interest  is  the  significantly  higher  prevalence  of  273 

A4396G  in  African  Americans.  African  American  women  274 

are  diagnosed  at  an  earlier  age  [28]  and  present  with  a  high-  275 

er  stage  at  diagnosis  [29].  Although  the  incidence  of  breast  276 

cancer  is  lower  [30],  except  for  very  young  women  [29],  277 

survival  also  is  lower  for  African  American  than  for  non-  278 

Hispanic  white  and  Hispanic  women  [31].  The  increased  279 

prevalence  of  the  A4396G  polymorphism,  particularly  if  280 

it  affects  IRF1  -mediated  tumor  suppression,  could  contrib-  281 

ute  to  these  observations  in  African  American  women,  but  282 
this  remains  to  be  established.  283 

284 

285 

286 
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Abstract 

Background:  The  poly  Q  polymorphism  in  AIBI  (amplified  in  breast  cancer)  gene  is  usually 
assessed  by  fragment  length  analysis  which  does  not  reveal  the  actual  sequence  variation.  The 
purpose  of  this  study  is  to  investigate  the  sequence  variation  of  poly  Q  encoding  region  in  breast 
cancer  cell  lines  at  single  molecule  level,  and  to  determine  if  the  sequence  variation  is  related  to 
AIBI  gene  amplification. 

Methods:  The  polymorphic  poly  Q  encoding  region  of  AIBI  gene  was  investigated  at  the  single 
molecule  level  by  PCR  cloning/sequencing.  The  amplification  of  A/8/  gene  in  various  breast  cancer 
cell  lines  were  studied  by  real-time  quantitative  PCR. 

Results:  Significant  amplifications  (5-23  folds)  of  AIBI  gene  were  found  in  2  out  of  9  (22%)  ER 
positive  cell  lines  (in  BT-474  and  MCF-7  but  not  in  BT-20,  ZR-75-1 ,  T47D,  BT483,  MDA-MB-361, 
MDA-MB-468  and  MDA-MB-330).  The  AIBI  gene  was  not  amplified  in  any  of  the  ER  negative  cell 
lines.  Different  passages  of  MCF-7  cell  lines  and  their  derivatives  maintained  the  feature  of  AIBI 
amplification.  When  the  cells  were  selected  for  hormone  independence  (LCCI)  and  resistance  to 
4-hydroxy  tamoxifen  (4-OH  TAM)  (LCC2  and  R27),  ICI  1 82,780  (LCC9)  or  4-OH  TAM,  KEO  and 
LY  I  17018  (LY-2),  AIBI  copy  number  decreased  but  still  remained  highly  amplified.  Sequencing 
analysis  of  poly  Q  encoding  region  of  AIBI  gene  did  not  reveal  specific  patterns  that  could  be 
correlated  with  AIBI  gene  amplification.  However,  about  72%  of  the  breast  cancer  cell  lines  had  at 
least  one  under  represented  (<20%)  extra  poly  Q  encoding  sequence  patterns  that  were  derived 
from  the  original  allele,  presumably  due  to  somatic  instability.  Although  all  MCF-7  cells  and  their 
variants  had  the  same  predominant  poly  Q  encoding  sequence  pattern  of 
(CAG)3CAA(CAG)9(CAACAG)3(CAACAGCAG)2CAA  of  the  original  cell  line,  a  number  of 
altered  poly  Q  encoding  sequences  were  found  in  the  derivatives  of  MCF-7  cell  lines. 

Conclusion:  These  data  suggest  that  poly  Q  encoding  region  of  AIBI  gene  is  somatic  unstable  in 
breast  cancer  cell  lines.  The  instability  and  the  sequence  characteristics,  however,  do  not  appear 
to  be  associated  with  the  level  of  the  gene  amplification. 
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Background 

While  predisposition  to  breast  cancer  is  largely  due  to 
mutations  in  high  penetrance  tumor  suppressor  genes 
such  as  BRCA1  and  BRCA2,  progression  of  cancer  is  the 
result  of  accumulation  of  genetic  alterations.  These  altera¬ 
tions  include  gene  amplifications,  microsatellite  instabili¬ 
ties,  loss  of  heterozygosity,  and  mutations  in  genes  that 
play  important  roles  in  signal  transduction  or  transcrip¬ 
tion  activation  pathways  leading  to  tumorigenesis.  Gene 
amplification  in  breast  cancer  was  found  in  several  chro¬ 
mosomal  locations  [1-4].  Among  them,  ErbB2  (or  HER-2/ 
neu)  amplification  strongly  correlates  with  steroid  recep¬ 
tor  negative  tumors  [5,6],  and  amplification  of 
AIB1  (amplified  in  breast  cancer  1)  gene  is  prevalent  in 
estrogen  receptor  (ER)  positive  tumors  [7,8].  The  AIB1 
gene  is  a  member  of  the  SRC-1  (steroid  receptor  coactiva¬ 
tor)  family  and  is  also  known  as  RAC3,  TRAM-1  or  ACTR 
[7,9,10].  It  is  located  at  chromosome  20ql2  region  and 
encodes  a  protein  of  1420  amino  acids  containing  bHLH- 
PAS  dimerization  domain,  a  hormone  receptor  interac¬ 
tion  domain,  a  CBP  interaction  domain,  and  histone 
acetyltransferase  domain  [11].  The  amplifications  and 
overexpression  of  AIB1  gene  were  found  to  be  a  common 
phenomenon  in  breast  cancer  cell  lines  and  primary 
breast  cancer  tissues  [12-15].  Since  AIB1  bridges  between 
nuclear  receptors  and  other  coactivators  or  the  transcrip¬ 
tional  machinery,  its  amplification  and  overexpression 
may  play  crucial  roles  in  the  development  of  breast  cancer 
and  may  potentially  have  influence  on  the  hormonal  pre¬ 
vention  and  treatment  for  breast  cancer. 

Toward  the  C-terminus  of  AIB1,  there  is  a  stretch  of  poly¬ 
glutamine  residues  that  are  encoded  by  polymorphic  CAG 
repeats.  The  expansion  of  CAG  repeats  in  poly  Q  contain¬ 
ing  proteins  underlies  a  number  of  neurodegenerative  dis¬ 
eases  [16,17].  Large  expansion  of  triplet  repeats  in  AIB1 
gene  does  not  occur,  presumably  due  to  the  frequent 
interruption  by  CAA  [18].  However  somatic  instability  by 
nucleotide  substitution  such  as  small  insertion  or  deletion 
does  occur  [18].  In  androgen  receptor  (AR),  the  length  of 
the  CAG  repeats  inversely  correlates  with  its  transcrip¬ 
tional  activity  [19,20].  Meanwhile  a  shorter  CAG  repeat  in 
AR  is  associated  with  a  higher  risk  of  an  aggressive  prostate 
cancer  phenotype  characterized  by  extraprostatic  exten¬ 
sion,  distant  metastases,  or  poor  histological  grade  [21]. 
In  the  case  of  AIB1,  it  is  not  clear  if  the  polymorphic 
length  of  poly  Q  affects  the  transactivation  activity  of 
AIB1.  AIB1  interacts  with  ER  in  a  ligand-dependent  man¬ 
ner  [7].  It  also  interacts  with  non-steroid  nuclear  receptors 
and  transcription  co-integrators  such  as  thyroid  and  retin¬ 
oid  receptors  and  CBP-dependent  transcription  com¬ 
plexes  [22,23].  Thus,  amplification  of  AIB1  gene  impacts 
on  both  estrogen  dependent  and  estrogen  independent 
mechanisms  leading  to  tumorigenesis  [24-26].  Although 
antiestrogens  are  the  most  common  type  of  endocrine 


therapy  in  breast  cancer  treatment,  acquired  resistance  can 
be  a  major  problem  in  clinical  management  of  initially 
responsive  breast  cancer  patients. 

Understanding  of  the  quantitative  and  qualitative  changes 
of  AIB1  gene  in  estrogen-independent  and  antiestrogen 
resistant  breast  cancer  cell  lines  may  help  in  the  selection 
of  steroid  or  non-steroid  antiestrogen  therapies.  Evalua¬ 
tion  of  AIB1  gene  amplification  in  previous  reports  is  per¬ 
formed  by  FISH  or  Southern  blot  analysis  [2,4,7].  In  this 
report,  we  use  the  real-time  quantitative  PCR  (Q-PCR) 
technique  to  assess  the  amplification  of  AIB1  gene  in  var¬ 
ious  breast  cancer  cell  lines  and  primary  breast  tumors. 
We  also  analyze  the  sequence  characteristics  and  instabil¬ 
ity  of  the  polymorphic  poly  Q  encoding  region  at  the  sin¬ 
gle  molecule  level  by  cloning  and  sequencing  of  the  DNA 
region  containing  CAG  repeats. 

Methods 

Samples  and  DNA  preparation 

Primary  breast  tumor  specimens  with  matching  normal 
breast  tissue  samples  were  obtained  from  Fu  Jen  Catholic 
University,  and  Cardinal  Tien  Hospital,  Taiwan,  after  sur¬ 
gical  removal  of  the  tumor  according  to  the  IRB  approved 
protocol.  The  ER  positive  breast  cancer  cell  lines  were 
obtained  from  Georgetown  University  Lombardi  Com¬ 
prehensive  Cancer  Center  Tissue  Culture  shared  resource 
and  American  Type  Cell  Culture.  A  total  of  25  cancer  and 
4  non-cancer  breast  tissue  cell  lines  were  studied.  MCF-7 
variants  include  different  passages,  MCF-7  pi 9,  MCF-7 
p72  and  MCF-7  derivatives:  LCC1  (selected  for  growth  in 
vitro  without  estrogens)  [27],  LCC2  (selected  from  LCC1 
by  treatment  with  non-steroid  antiestrogen  4-OH  TAM) 
[28],  LCC9  (selected  from  LCC1  by  treatment  with  steroid 
antiestrogen  ICI  182,780)  [29],  LY-2  (resistant  to  4-OH 
TAM,  KEO  and  LY  117018)  and  R27  (resistant  to  4-OH 
TAM).  AK-47  is  derived  from  parental  ER  positive  cell  line 
ZR-75-1  with  the  loss  of  expression  of  ER.  LCC6  is  a  more 
aggressive  form  of  MDA-MB-435  [30].  A1N4  is  a  normal 
breast  cell  line  that  is  ER  negative.  DNA  from  tissues, 
blood  and  cell  cultures  was  extracted  by  salting  out 
method  [31]. 

Preparation  of  standard  DNA  for  quantitative  PCR 

A  region  of  439  bp  from  exon  5  of  AIB1  gene  was  ampli¬ 
fied  with  the  forward  primer;  5'-CAAGCGATCAAATGAG- 
GGTAG-3'  and  the  reverse  primer;  5'- 
CATTGTrTCATATCTCTGGCG-3 ' .  A  fragment  of  85  bp 
from  3'  untranslated  region  of  j32-micro globulin  gene  (/?2- 
M)  was  amplified  with  the  forward  primer;  5'-TGCT- 
GTCTCCATGTlTGATGTATCT-3'  and  the  reverse  primer; 
5 '-TCTCTGCTCCCCACCTCTAAGT-3 '  [32-34].  These  PCR 
products  were  cloned  into  the  pCR  2.1-TOPO  vector  (Inv- 
itrogen).  The  plasmid  DNA  was  isolated  and  quantified 
using  the  DU640  Spectrophotometer  (Beckman,  Fuller- 
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ton,  CA,  USA).  The  copy  numbers  were  calculated  from 
absorbance  at  260  nm  and  based  on  the  molecular  weight 
of  the  resulting  plasmid.  The  plasmid  DNAs  were  serially 
diluted  over  four  logs  to  establish  the  standard  curve  giv¬ 
ing  a  range  from  400,000  to  40  copies/pl.  In  additional  set 
of  experiments  the  standard  curve  was  constructed  using 
genomic  DNA  prepared  as  a  pool  of  equal  amounts  of 
blood  DNA  from  7  control  individuals  with  normal  AIB1 
copy  number.  'Normal'  genomic  DNA  (100  ng/pl)  was 
diluted  in  water  over  four  logs.  Since  1  ng  of  genomic 
DNA  contains  approximately  330  copies  of  a  single  copy 
gene,  five  standards  used  range  from  33000  to  3.3  copies/ 

hi- 

Real-time  quantitative  PCR  (RT  Q-PCR) 

In  real-time  Q-PCR  analysis,  the  primers  used  were  5'- 
GAGTTTCCTGGACAAATGAG-3 '  (forward)  and  5'-CATT- 
GTTTCATATCTCTGGCG-3 '  (reverse)  for  AIB1  gene  (Exon 
5),  and  the  same  primers  as  used  for  standard  DNA  prep¬ 
aration  for  JB2-M  gene,  yielding  134  bp  and  85  bp  PCR 
products,  respectively.  The  TaqMan  probes  were  FAM-5' 
GCCGTATGTTGATGAAAACACCACA  3'-TAMRA  and  VIC- 
5'  TTGCTCCACAGGTAGCTCTAGGAGG  3'-TAMRA,  for 
AIB1  and  J32-M  gene  respectively,  each  labeled  with  FAM 
or  VIC  (reporter  dye)  at  the  5'  end  and  TAMRA  (quencher 
dye)  at  the  3'  end.  Each  10  pi  real  time  Q-PCR  reaction 
mixture  contained  1  x  TaqMan  Universal  PCR  Master  Mix 
(Applied  Biosystems,  Foster  City,  CA),  10  ng  of  genomic 
DNA,  0.3  pM  of  each  primer,  and  0.1  pM  probe.  The  actin 
gene  was  also  used  as  a  reference.  However,  since  the  actin 
gene  has  multiple  homologous  copies,  the  data  presented 
here  were  referenced  to  J32-M  gene.  The  amplification  was 
carried  out  according  to  the  conditions  suggested  by  the 
manufacturer  (initial  denaturation  at  95  °C  for  10  min 
and  40  cycles  of  95°  C  for  15  s  and  60  °C  for  1  min)  using 
an  ABI  Prism  7700  Sequence  Detection  System  (Applied 
Biosystems,  Foster  City,  CA).  Each  measurement  was  per¬ 
formed  in  triplicate  and  the  threshold  cycle  numbers  ( C T) 
were  measured.  The  copy  number  was  generated  from  the 
Cr  value  and  standard  curve  according  to  previously 
described  procedures  [32-34]. 

Cloning  and  sequencing 

The  poly  Q  containing  fragment  was  amplified  by  the  for¬ 
ward  primer  F:  5'  GTCTTATACCTGGTGTATTG  3'  and  the 
reverse  primer  R:  5'  CTGGGGGAAGCAGTCACATTAG  3', 
yielding  a  PCR  product  of  3 14  bp.  The  high  fidelity  ampli¬ 
fication  was  carried  out  in  a  30  pi  reaction  mixture  con¬ 
taining  10  ng  of  genomic  DNA,  0.2  pM  of  each  primer,  1 
x  HF  2  PCR  buffer,  dNTPs,  and  Advantage-HF  2  polymer¬ 
ase  according  to  the  manufacturer's  recommendation 
(Clontech  Faboratories,  Palo  Alto,  CA).  After  1  min  of  ini¬ 
tial  denaturation  at  94 °C,  the  DNA  was  amplified  by  30 
cycles  of  45  s  at  95  ° C,  45  s  at  55  °  C  and  45  s  at  72 ° C,  fol¬ 
lowed  by  a  final  extension  at  72  °C  for  5  min.  The  PCR 


products  were  purified  and  cloned  into  pCR2.1-TOPO 
(Invitrogen)  vector  according  to  the  manufacturer's  proto¬ 
col.  At  least  8  clones  from  each  sample  were  picked  for 
sequencing  using  BigDye  sequencing  kit  and  analyzed  on 
an  ABI  377  DNA  Sequencer  (Applied  Biosystems,  Foster 
City,  CA).  Two  primers,  F  and  F2  (5'  AGCAGGGl  11 1CT- 
TAATGCTC  3')  were  used  for  sequencing  and  loading  of 
reactions  onto  alternate  lanes  for  easy  tracking.  The 
sequence  results  were  analyzed  using  sequence  analysis 
software  version  3.4. 

Results 

Amplification  ofAIBI  gene 

Real-time  Q-PCR  analysis  allows  the  measurement  of 
actual  copy  number  of  AIB1  gene  using  a  single  copy  /?2- 
microglobulin  gene  as  a  reference.  From  the  threshold  cycle 
number  and  the  standard  curve,  the  ratio  of  the  copy 
number  of  AIB1  gene  to  that  of  J32-M  gene  can  be  calcu¬ 
lated.  This  ratio  can  be  used  as  a  measure  of  the  amplifica¬ 
tion  of  the  AIB1  gene.  The  average  copy  number  ratio  of 
the  AIBl//32-M  in  the  blood  samples  from  48  age  matched 
control  individuals  is  determined  to  be  1.16  ±  0.38.  An 
AIB1/J32-M  ratio  above  2  SD  of  the  mean  (1.16  +  2  x  0.38 
=  1.92)  is  defined  as  truly  amplified.  In  addition,  all  meas¬ 
urements  were  repeated  using  a  pool  of  normal  genomic 
DNAs  for  standard  curve  construction.  The  results 
obtained  using  both  methods  were  practically  identical. 

We  first  evaluate  the  amplification  of  AIB1  gene  in  26  pri¬ 
mary  breast  tumors  (13  ER  positive  and  13  ER  negative) 
and  corresponding  surrounding  normal  breast  tissue  sam¬ 
ples.  AIB1  gene  was  found  to  be  amplified  in  1  ER  positive 
tumor  sample  that  constitutes  3.8%  of  total  or  7.6%  of  ER 
positive  tumors.  This  result  is  consistent  with  previous 
report  [3,7,13]. 

As  shown  in  Table  1,  9  out  of  29  cell  lines  had  elevated 
AIB1  at  2SD  above  the  mean.  All  of  them  were  ER  positive. 
AIB1  gene  was  not  amplified  in  7  ER  positive  cell  lines: 
BT20,  ZR-75-1,  T47D,  BT483,  MDA-MB-361,  MDA-MB- 
468  and  MDA-MB-330.  None  of  the  13  ER  negative  cell 
lines  showed  significant  AIB1  amplification.  Some  cell 
lines  are  derivatives  of  others.  For  example,  AK-47  is 
derived  from  ER  positive  ZR-75-1  cell  line.  In  AK-47  cells, 
loss  of  ER  expression  did  not  have  any  effect  on  AIB1  copy 
number.  The  amplification  of  AIB1  gene  in  different  pas¬ 
sages  of  ER  positive  MCF-7  cell  lines  remains  at  high  levels 
of  18-23  fold  of  control.  Loss  of  estrogen  dependence  in 
LCC1  is  accompanied  by  moderate  decrease  in  AIB1  gene 
amplification  (13.5  fold  in  LCC1  versus  22.2  fold  in  MCF- 
7).  The  level  of  AIB1  gene  amplification  was  reduced  to 
14.6  fold  when  the  estrogen-independent  cells  became 
resistant  to  antiestrogen  4-OH  TAM  treatment  (LCC2). 
Similar  decrease  in  amplification  (15.6  fold  of  control)  of 
AIB1  gene  was  observed  in  estrogen-independent  cell  line 
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Table  I:  The  ER  status,  AIBI  amplification,  and  poly  Q  encoding 
sequences  in  breast  cell  lines 


Cell  line  0 

ER  status 

AIBI  copy  number 
ratio 

(tumor/normal) 

Number  of  s 
equence 
patterns 

MCF-IOA 

_ 

l.l 

3 

MCF-IOA  neo 

- 

1.3 

3 

Al  N4b 

- 

2.0 

2 

AK-47 

- 

1.4 

3 

MDA-MB435 

- 

1.3 

4 

LCC6 

- 

2.0 

7 

MDA-MBI 57 

- 

l.l 

6 

MDA-MBI 34V 

- 

1.4 

2 

MDA-MB23 1  N 

- 

1.2 

2 

HBLI00 

- 

1.2 

1 

ZR-75-30 

- 

1.7 

1 

HS  578T 

- 

1.9 

2 

HS  578BST 

- 

l.l 

2 

MCF7 

+ 

22.2 

1 

MCF-7  PI 9 

+ 

18.7 

3 

MCF-7  P72 

+ 

22.8 

3 

LY-2 

+ 

19.9 

2 

R27 

+ 

19.4 

2 

LCCI 

+ 

13.5 

3 

LCC2 

+ 

14.6 

2 

LCC9 

+ 

15.6 

2 

BT474 

+ 

4.9 

2 

BT-483 

+ 

1.9 

4 

BT20 

+ 

1.7 

2 

MDA-MB468 

+ 

1.7 

2 

ZR-75- 1 

+ 

1.5 

2 

MDA-MB36I 

+ 

1.3 

2 

T47D 

+ 

1.0 

5 

MDA-MB330 

+ 

0.9 

2 

°LCCI :  estrogen  independent  and  responsive  which  is  selected  for 
growth  in  vivo  without  estrogens;  LCC2:  selected  from  LCCI  by 
treatment  with  4-OH  TAM;  LCC9:  selected  from  LCCI  by  treatment 
with  ICI  182780,  resistant  to  ICII82780  and  4-OH  TAM:,  LY-2: 
selected  for  resistance  to  4-OH  TAM,  KEO  and  LY  I  17018;  R27:  able 
to  grow  in  the  presence  of  4-OH  TAM.  LCC6  is  the  more  aggressive 
variant  of  MD-MB435.  AK-47  is  the  variant  derived  from  ZR-75- 1 . 

A I  N4  is  a  normal  breast  cell  line. 
bCell  lines  listed  in  Table  2  are  in  bold. 

that  has  gained  antiestrogen  resistance  to  both  non-ster¬ 
oid,  4-OH  TAM,  and  steroid,  ICI  182,780  antiestrogens 
(LCC9)  (Table  1). 

Somatic  instability  of  poly  Q  encoding  region  of  AIBI  gene 
in  breast  cancer  cell  lines 

The  polymorphic  poly  Q  encoding  region  of  AIBI  con¬ 
tains  CAG  repeat  that  is  frequently  interrupted  by  CAA's. 
The  poly  Q  region  is  part  of  the  histone  acetyltransferase 
domain.  It  is  also  where  the  recruitment  and  interactions 
with  other  components  of  the  transcription  activator  com¬ 
plex  takes  place.  In  order  to  investigate  if  qualitative  alter¬ 
ation  in  this  region  accompanied  the  quantitative  change 
of  AIBI  gene  in  breast  cancer  cell  lines,  we  cloned  and 
sequenced  the  poly  Q  encoding  region  of  the  gene.  The 


cloning/sequencing  technique  resolved  the  heterogene¬ 
ous  poly  Q  encoding  sequences  into  distinct  sequences, 
thus  allowing  the  analysis  at  the  single  molecule  level.  At 
least  8  clones  from  each  cell  line  were  selected  and 
sequenced.  Theoretically,  there  should  be  only  2  distinct 
sequence  patterns  if  the  cell  line  is  heterozygous  for  AIBI 
allele  and  one  distinct  sequence  pattern  if  it  is 
homozygous.  However,  18/25  (72%)  (data  partially 
shown  in  Table  2)  of  the  cell  lines  contain  at  least  one 
poly  Q  encoding  sequence  pattern  that  represents  less 
than  20%  of  the  sequenced  clones  of  the  cell  line.  These 
results  suggest  that  the  under-represented  sequences  prob¬ 
ably  arise  from  the  parental  sequence  by  somatic  muta¬ 
tion.  Indeed  these  rare  sequences  differ  from  their 
parental  sequence  by  one  base  pair  substitution  (CAG  to 
CAA)  or  by  insertion  or  deletion  of  CAGs.  The  high  degree 
(72%  of  the  cell  lines)  of  somatic  instability  is  probably 
characteristic  for  cancer  cell  lines  since  it  only  occurs  in 
less  than  5%  (2/43)  of  the  normal  controls.  We  analyzed 
normal  A1N4  cell  line  at  two  different  times.  The  first 
time,  ten  clones  of  A1N4  cell  line  were  sequenced.  Pattern 
2  (Table  2),  (CAG)6CAA(CAG)9(CAACAG)3(CAACAG- 
CAG)2CAA,  was  found  in  4  clones,  and  pattern  17, 
(CAG)4CAA(CAG)9(CAACAG)3(CAACAGCAG)2CAA  was 
found  in  6  clones.  The  second  time,  31  clones  were 
sequenced.  Sixteen  had  pattern  2  and  15  had  pattern  17. 
There  was  no  occurrence  of  "extra"  poly  Q  encoding 
sequence.  These  results  suggest  that  the  occurrence  of  rare 
sequences  is  not  due  to  cloning/PCR  artifact. 

Association  between  poly  Q  length  or  its  specific  encoding 
sequence  with  AIBI  amplification  was  not  recognized 
(Table  2).  Somatic  instability  occurred  in  all  variants  of 
MCF-7  cell  line,  although  they  all  maintain  the  parental 
allele  as  the  predominant  coding  sequences  (pattern  1). 
Two  new  sequences  arise  by  insertion  of  2  and  3  CAGs  in 
passage  19.  Another  two  new  sequences  occur  in  passage 
72  by  deletion  of  1  and  2  CAG  repeats.  Similar  somatic 
mutations  occur  in  cell  lines  LCCI,  LCC2,  LCC9,  LY-2  and 
R27.  These  mutations  seem  to  occur  randomly  and  inde¬ 
pendently  in  each  cell  line.  There  is  not  any  single 
sequence  that  occurs  more  frequently  than  others,  except 
pattern  5,  which  occurs  3  times  by  losing  1  CAG  directly 
from  the  parental  sequence. 

AK47  was  derived  from  ZR-75- 1  by  losing  its  ER  activity. 
During  the  establishment  of  the  cell  line,  additional 
somatic  mutations  occurred  in  the  polyglutamine  region 
(patterns  11  and  16).  The  poly  Q  encoding  sequence  of 
AIBI  gene  seems  to  be  quite  unstable  in  MDA-MB435  cell 
line.  It  has  4  distinct  poly  Q  encoding  sequence  patterns 
with  pattern  9  as  the  predominant  one.  Its  variant  LCC6 
had  a  total  of  7  different  sequence  patterns.  These  data  are 
consistent  with  the  genomic  instability  that  is  characteris¬ 
tic  for  cancer  cells.  Although  poly  Q  encoding  sequence 
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patterns  do  not  seem  to  directly  link  to  AIB 1  gene  amplifi¬ 
cation,  it  is  possible  that  the  alteration  in  poly  Q  length 
affects  protein-protein  interaction,  thus,  the  transactiva¬ 
tion  activity  of  AIB  1.  While  most  alterations  do  not  change 
poly  Q  length  significantly,  rare  sequence  pattern  in  LCC2 
with  much  shorter  (only  14  repeats)  poly  glutamine  tract 
may  affect  the  co-transactivating  activity  of  AIB1  gene. 

Discussion 

Genetic  and  clinical  phenotypic  heterogeneity  is  the 
prominent  characteristic  of  breast  cancer.  Multiple  genetic 
alterations  contribute  to  breast  cancer  development  and 
progression  [35,36].  The  occurrence  of  DNA  amplifica¬ 
tions  in  breast  cancer  had  been  studied  by  Southern  blot 


[1],  FISH  (fluorescence  in  situ  hybridization)  [4]  and 
CGH  methods  (comparative  genomic  hybridization)  [37- 
39].  We  developed  real  time  quantitative  PCR  method  to 
more  accurately  assess  the  amplification  of  AIB1  gene  in 
breast  cancer  cell  lines.  Amplification  of  AIB1  in  breast 
cancer  cell  lines;  BT-474  and  MCF-7  were  first  reported  by 
Guan  et  al.  [3].  By  FISH  analysis,  Anzick  et  al.  [7]  observed 
>20  fold  amplification  of  AIB1  gene  in  three  ER  positive 
breast  cancer  cell  lines  (BT-474,  MCF-7  and  ZR-75-1)  and, 
to  a  lesser  extent,  in  10%  primary  breast  tumors.  In  this 
study  we  did  not  detect  significant  amplification  in  ZR-75- 
1  cell  line.  The  discrepancy  may  be  explained  by  a  differ¬ 
ent  source  of  the  cell  line  or  by  spontaneous  change  of  the 
cell  line  during  passages.  In  addition,  FISH  analysis  is 


Table  2:  Poly  Q  sequence  patterns  and  AIB  I  amplification  level  in  MCF7  and  its  variants 


Cell  line0 

amplification 

Sequence  patterns 

Pattern 

(Q)n 

Frequency 

MCF-7 

22.3 

(CAG)3CAA(CAG)9(CAACAG)3(CAACAGCAG)2CAA 

1 

26 

7/7 

MCF-7  PI 9 

18.7 

(CAG)3CAA(CAG)9(CAACAG)3(CAACAGCAG)2CAA 

1 

26 

7/9 

(CAG)6CAA(CAG)9(CAACAG)3(CAACAGCAG)2CAA 

2 

29 

1/9 

(CAG)sCAA(CAG)9(CAACAG)3(CAACAGCAG)2CAA 

3 

28 

1/9 

MCF-7  P72 

22.8 

(CAG)3CAA(CAG)9(CAACAG)3(CAACAGCAG)2CAA 

1 

26 

5/7 

(CAG)3CAA(CAG)8(CAACAG)3(CAACAGCAG)2CAA 

4 

25 

1/7 

(CAG)3CAA(CAG)7(CAACAG)3(CAACAGCAG)2CAA 

5 

24 

1/7 

LCC  1 

13.5 

(CAG)3CAA(CAG)9(CAACAG)3(CAACAGCAG)2CAA 

1 

26 

7/9 

(CAG)3CAA(CAG)8(CAACAG)3(CAACAGCAG)2CAA 

4 

25 

1/9 

(CAG)3CAA(CAG)I4(CAACAG)2(CAACAGCAG)2CAA 

6 

29 

1/9 

LCC2 

14.6 

(CAG)3CAA(CAG)9(CAACAG)3(CAACAGCAG)2CAA 

1 

26 

4/5 

CAG(CAACAG)3(CAACAGCAG)2CAA 

7 

14 

1/5 

LCC9 

15.6 

(CAG)3CAA(CAG)9(CAACAG)3(CAACAGCAG)2CAA 

1 

26 

7/8 

(CAG)3CAA(CAG)2CAA(CAG)9(CAACAG)3(CAACAGCAG)2CAA 

8 

29 

1/8 

LY-2 

19.9 

(CAG)3CAA(CAG)9(CAACAG)3(CAACAGCAG)2CAA 

1 

26 

9/10 

(CAG)6CAA(CAG)8(CAACAG)3(CAACAGCAG)2CAA 

9 

28 

1/10 

R27 

19.4 

(CAG)3CAA(CAG)9(CAACAG)3(CAACAGCAG)2CAA 

1 

26 

8/9 

(CAG)3CAA(CAG)8(CAACAG)3(CAACAGCAG)2CAA 

4 

25 

1/9 

MDA-MB435 

1.3 

(CAG)3CAA(CAG)2CAA(CAG) , ,  (CAACAG)2(CAACAGCAG)2CAA 

10 

29 

2/8 

(CAG)6CAA(CAG)8(CAACAG)3(CAACAGCAG)2CAA 

9 

28 

3/8 

(CAG)3CAA(CAG)9(CAACAG)3(CAACAGCAG)2CAA 

1 

26 

2/8 

(CAG)6CAA(CAG)8(CAACAG)2(CAACAGCAG)2CAA 

1 1 

26 

1/8 

LCC6 

2.0 

(CAG)3CAA(CAG)2CAA(CAG) ,  0(CAACAG)2(CAACAGCAG)2CAA 

12 

28 

1/9 

(CAG)3CAA(CAG)2CAA(CAG)6(CAACAG)3(CAACAGCAG)2CAA 

13 

26 

1/9 

(CAG)6CAA(CAG)8(CAACAG)3(CAACAGCAG)2CAA 

9 

28 

3/9 

(CAG)3CAA(CAG)2CAA(CAG)io(CAACAG)3(CAACAGCAG)2CAA 

14 

30 

1/9 

(CAG)6CAA(CAG)8(CAACAG)3(CAACAGCAG)2CAA 

15 

28 

1/9 

(CAG)3CAA(CAG)2CAA(CAG) , ,  (CAACAG)2(CAACAGCAG)2CAA 

10 

29 

1/9 

(CAG)3CAA(CAG)2CAACAGCAA(CAG)8(CAACAG)2(CAACAGCAG)2CAA 

16 

28 

1/9 

ZR-75-1 

1.5 

(CAG)6CAA(CAG)9(CAACAG)3(CAACAGCAG)2CAA 

2 

29 

6/7 

(CAG)sCAA(CAG)9(CAACAG)3(CAACAGCAG)2CAA 

3 

28 

1/7 

AK-47 

1.4 

(CAG)6CAA(CAG)9(CAACAG)3(CAACAGCAG)2CAA 

2 

29 

2/8 

(CAG)6CAA(CAG)8(CAACAG)3(CAACAGCAG)2CAA 

1 1 

28 

5/8 

(CAG)sCAA(CAG)8(CAACAG)3(CAACAGCAG)2CAA 

15 

27 

1/8 

AIN4 

2.0 

(CAG)6CAA(CAG)9(CAACAG)3(CAACAGCAG)2CAA 

2 

29 

20/41 

(CAG)4CAA(CAG)9(CAACAG)3(CAACAGCAG)2CAA 

17 

27 

21/41 

°LCC  I :  estrogen  independent  and  responsive  which  is  selected  for  growth  in  vivo  without  estrogens;  LCC2:  selected  from  LCC I  by  treatment  with 
4-OH  TAM;  LCC9:  selected  from  LCC  I  by  treatment  with  ICI  182780,  resistant  to  ICI 182780  and  4-OH  TAM:,  LY-2:  selected  for  resistance  to  4- 
OH  TAM,  KEO  and  LY  I  17018;  R27:  able  to  grow  in  the  presence  of  4-OH  TAM.  LCC6  is  the  more  aggressive  variant  of  MD-MB435.  AK-47  is  the 
variant  derived  from  ZR-75-1.  AIN4  is  a  normal  breast  cell  line.  The  parental  sequence  patterns  are  in  bold. 
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restricted  to  a  few  cells,  whereas  real  time  qPCR  analysis 
measures  the  gene  in  the  overall  DNA  extract.  Glaeser  et 
al.  [2]  used  the  quantitative  differential  PCRto  determine 
the  amplification  level  of  AIB1  and  found  no  amplifica¬ 
tion  in  breast  or  endometrial  carcinomas.  These  methods 
did  not  give  actual  copy  numbers  of  the  gene.  In  this 
study,  we  used  real-time  Q-PCR  to  determine  the  level  of 
AIB1  amplification.  The  ability  of  real-time  Q-PCR  to 
detect  the  fluorescent  signal  from  degraded  sequence  spe¬ 
cific  TaqMan  probe  at  the  very  beginning  period  of  expo¬ 
nential  stage  offered  an  accurate  way  of  DNA 
quantification.  When  compared  to  CGH,  FISH  and  South¬ 
ern  blot  analysis,  this  method  has  the  advantage  of  high 
sensitivity,  reproducibility,  and  efficiency. 

If  only  the  original  cell  lines  were  counted,  the  AIB1  gene 
was  amplified  in  2  out  of  9  ER  positive  and  none  of  10  ER 
negative  cell  lines.  Higher  degree  (22%)  of  AIB1  amplifi¬ 
cation  in  ER  positive  breast  cancer  cell  lines  may  suggest 
the  association  between  AIB 1  gene  amplification  and  ER 
status.  This  is  further  supported  by  the  observation  that 
the  AIB1  amplification  moderately  decreased  when  cells 
became  ER  independent  as  LCC1,  LCC2  and  LCC9  con¬ 
sistent  with  the  role  of  AIB1  in  ER-dependent  signaling. 
All  MCF-7  variant  cell  lines  maintain  high  level  of  AIB1 
gene  amplification  of  the  parental  cells.  Additional  gain  of 
resistance  to  an  antiestrogen,  ICI  182,780,  does  not  have 
significant  effect  on  AIB1  amplification  (from  LCC1  to 
LCC9).  Similarly,  resistance  to  4-OH  TAM  is  not  consist¬ 
ently  accompanied  with  the  change  in  AIB1  gene  amplifi¬ 
cation.  LY-2,  R27,  and  LCC2  were  all  selected  from 
estrogen  dependent  MCF-7  cells  against  high  dose  of  4- 
OH  TAM.  AIB1  gene  amplification  in  LY-2  and  R27  cell 
lines  remained  almost  unchanged,  whereas  in  LCC2,  AIB1 
gene  amplification  is  moderately  decreased.  These  data 
suggest  that  resistance  to  4-OH  TAM  does  not  necessarily 
affect  AIB  1  gene  amplification.  It  should  be  noted  that  in 
LCC2  there  is  a  short  poly  Q  containing  mutant  AIB  1  and 
lower  expression  of  the  gene  may  be  compensated  by  the 
increase  in  co-transactivation  activity  of  the  mutant  pro¬ 
tein.  Our  observation  of  variations  in  AIB1  gene  amplifi¬ 
cation  in  various  derivatives  of  MCF-7  cell  dines  is 
consistent  with  the  previous  report  in  which  several  MCF- 
7  sublines  were  shown  to  have  the  capacity  to  generate 
clonal  heterogeneity.  This  represents  an  important  selec¬ 
tive  advantage  in  MCF-7  in  leading  to  aggressive  and  met¬ 
astatic  forms  of  the  disease  [40]. 

Besides  the  quantitative  regulation  of  AIB  1  gene  in  breast 
cancer  cell  lines,  the  AIB1  gene  contains  CAG  repeat 
region  which  is  a  target  for  genetic  instability  in  tumor 
progression.  Large  expansion  of  triplet  repeat  occurs  in 
various  neuro degenerative  diseases  [41-43].  These  abnor¬ 
mal  proteins  form  large  aggregates  that  have  been  shown 
to  tie  up  transcription  factors  that  bind  poly  Q  such  as 


CREB  [44].  The  poly  Q  tract  in  the  androgen  receptor  (AR) 
gene  is  unique  in  that  the  large  expansion  of  poly  Q 
encoding  CAG  repeat  causes  the  X-linked  spinal  bulbar 
muscular  atrophy  (SBMA,  or  Kennedy  Disease)  [19,20], 
but  short  poly  Q  of  AR  is  correlated  with  hormone- 
dependent  transactivation  [19,20]  and  more  aggressive 
form  of  cancer  [21,45].  AIB1  shares  several  structural/ 
functional  similarities  with  AR.  Both  genes  are  involved  in 
nuclear  receptor  mediated  regulation  of  gene  expression. 
Due  to  frequent  interruption  with  CAA's,  large  expansion 
of  the  triplet  was  not  observed  in  AIB1  gene.  This  region, 
however,  was  quite  unstable  as  evidenced  by  frequent 
CAA/CAG  changes  and  small  insertions  and  deletions. 
The  PCR/ cloning  strategy  allows  us  to  investigate  the  pol¬ 
ymorphic  poly  Q  encoding  region  at  single  molecule 
level.  Since  we  only  sequenced  a  small  number  of  clones 
from  each  individual,  we  cannot  exclude  the  possibility 
that  the  under-represented  alleles  may  be  lost  in  PCR/ 
cloning/sequencing  process.  Several  distinct  poly  Q 
encoding  sequence  patterns  were  observed  in  LCC-6, 
T47D,  and  MDA-MB157.  T47D  is  a  cell  line  of  notable 
genetic  instability  that  was  observed  in  the  estrogen  recep¬ 
tor  gene  [46].  LCC6  cell  line  which  formed  ascites  was  a 
more  aggressive  form  of  MDA-MB435  [30].  Since  various 
rare  poly  Q  encoding  sequences  seem  to  arise  randomly 
and  independently  in  different  cell  lines  regardless  of  the 
AIB1  gene  amplification  levels,  we  attribute  these  somatic 
mutations  to  genetic  mutability  in  cancer  cells. 

Conclusion 

The  poly  Q  encoding  sequence  of  AIB  1  gene  is  genetically 
unstable  and  is  an  easy  target  for  somatic  mutations  in 
cancer  cells.A7Bl  gene  amplification  occurs  in  only  a 
small  fraction  of  ER  positive  primary  breast  tumors  and 
breast  cancer  cell  lines.  AIB1  gene  amplification  has  not 
been  found  in  ER  negative  primary  tumor  or  breast  cancer 
cell  lines.  Gain  of  estrogen  independence  and  resistance  to 
steroid  antiestrogen  may  be  accompanied  by  moderate 
decrease  of  AIB  1  gene  amplification. 

Abbreviations 

4-OH  TAM  =  4-hydroxy  tamoxifen;  ACTR  =  Activator  of 
Retinoid  and  Thyroid  Receptors;  AIB1  =  Amplified  in  Breast 
Cancer  gene  1;  bHLH-PAS  =  basic  helix-loop-helix  -Per- 
ARNT-Sim;  J32-M  =  fi2-micro globulin  gene;  CBP  =  CREB 
binding  protein;  CREB  =  cAMP-responsive  element  bind¬ 
ing  protein;  ER  =  estrogen  receptor;  FAM  =  6-Carboxy-Flu- 
orescein;  RAC3  =  Receptor- Associated  Co-activator  3; 
TRAM-1  =  Thyroid  Hormone  Receptor  Activator 
Moleculel;  TAMRA  =  6-Carboxy-Tetramethyl  rhodamine 
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Abstract 

Hormone-dependent  estrogen  receptor  (ER)-positive  breast  cancer  cells  may  adapt  to  low 
estrogen  environments  such  as  produced  by  aromatase  inhibitors.  In  many  instances,  cells 
become  insensitive  to  the  effects  of  estrogen  but  may  still  retain  dependence  on  ER.  We  have 
investigated  the  expression,  function,  and  activation  of  ERa  in  two  endocrine-resistant  MCF-7 
models  to  identify  mechanisms  that  could  contribute  to  resistance.  While  MCF-7/LCC1  cells  are 
partially  estrogen  dependent,  MCF-7/LCC9  cells  are  fully  estrogen  insensitive  and  fulvestrant  and 
tamoxifen  resistant.  In  both  MCF-7/LCC1  and  MCF-7/LCC9  cell  lines,  high  expression  of  ERa  was 
associated  with  enhanced  binding  to  the  trefoil  factor  1  (TFF1)  promoter  in  the  absence  of 
estrogen  and  increased  transcription  of  TFF1  and  progesterone  receptor.  In  contrast  to  the 
observations  derived  from  hypersensitive  and  supersensitive  models,  these  cells  were  truly 
estrogen  independent;  nevertheless,  removal  of  ERa  by  siRNA,  or  fulvestrant,  a  specific  ER 
downregulator,  inhibited  growth  indicating  dependence  on  ERa.  In  the  absence  of  estrogen, 
neither  ERa  Ser118  nor  Ser167  were  phosphorylated  as  frequently  found  in  other  ligand- 
independent  cell  line  models.  Addition  of  estrogen  activated  ERa  Ser118  in  MCF-7  and  LCC1  cells 
but  not  in  LCC9  cells.  We  suggest  that  the  estrogen-independent  growth  within  these  cell  lines  is 
accounted  for  by  high  levels  of  ERa  expression  driving  transcription  and  full  estrogen 
independence  explained  by  lack  of  ERa  activation  through  Ser118. 
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Introduction 

Estrogen  receptor  a  (ERa)  is  a  major  growth  regulator 
for  many  breast  cancers  and  has  provided  an 
exploitable  target  for  therapy  (Ali  &  Coombes  2002). 
Estrogen  binding  to  ERa  promotes  conformational 
changes  in  the  receptor  leading  to  dimerization  and 
attachment  to  DNA,  generally  at  the  site  of  conserved 
estrogen  response  elements  in  the  promoter  regions  of 
target  genes  (Ali  &  Coombes  2002).  Functional 
regulation  of  ERa  is  additionally  mediated  via 
phosphorylation  of  key  residues  in  the  activation 
function  1  (AF-1)  domain  of  ERa  including  Ser118 
and  Ser167  and  these  influence  both  DNA  binding  and 


recruitment  of  cofactor  molecules  (reviewed  in 
Lannigan  2003).  The  activation  of  ER  involves 
crosstalk  with  other  growth  factor- signaling  pathways. 
There  is  extensive  evidence  that  activation  of  the 
mitogen-activated  protein  kinase  (MAPK)- signaling 
cascade  and  the  phosphoinositol  3  kinase  (PI3-K) 
pathway  phosphorylate  ERa  at  Ser118  and  Ser167,  via 
extracellular  signal-regulated  kinase  (ERK)l/2  and  Akt 
respectively  (Bunone  et  al.  1996,  Martin  et  al.  2000, 
Lannigan  2003).  Transcriptional  activation  of  ERa 
then  involves  a  dynamic  process  where  large  transcrip¬ 
tion  complexes  incorporating  co-activator  proteins  are 
assembled  in  an  ordered  and  combinatorial  manner 
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(Glass  &  Rosenfeld  2000,  Metivier  et  al.  2003). 
Well-defined  estrogen-regulated  genes  include  trefoil 
factor  1  (TFFl)/pS2  (Masiakowski  et  al.  1982, 
Jakowlew  et  al.  1984)  and  progesterone  receptor 
(PGR;  Nardulli  et  al.  1988). 

While  tamoxifen  has  been  the  established  form  of 
treatment  for  ER-positive  breast  cancers  for  more  than 
20  years,  other  anti-estrogen  strategies,  notably 
aromatase  inhibitors  (Johnston  &  Dowsett  2003)  and 
selective  estrogen  downregulators  (SERDs),  are 
increasingly  being  used  (Robertson  2002).  Despite 
initial  responsiveness  to  these  agents,  most  tumors 
eventually  recur  with  acquired  resistance  (Clarke  et  al. 

2001,  2003).  Multiple  mechanisms,  dependent  on  the 
form  of  endocrine  treatment,  are  involved  in  the 
development  of  resistance  and,  in  many  cases,  these 
mechanisms  remain  unclear.  During  the  acquisition  of 
endocrine  resistance,  progressive  changes  are  fre¬ 
quently  observed,  with  ER-positive  breast  cancer 
cells  progressing  in  a  stepwise  manner  from  a  fully 
estrogen-sensitive  phenotype  to  an  estrogen- sensitive, 
but  no  longer  dependent  phenotype,  to  a  fully  resistant 
phenotype  (Clarke  et  al.  2001,  2003). 

With  the  increasing  clinical  use  of  aromatase 
inhibitors,  such  as  letrozole,  anastrazole,  and 
exemestane  which  act  by  inhibiting  estrogen 
synthesis  (Johnston  &  Dowsett  2003),  there  has 
been  great  interest  in  how  breast  cancer  cells  can 
adapt  to  low  estrogen  environments  and  become 
resistant  to  the  effects  of  these  drugs.  In  most  cases 
of  acquired  anti-estrogen  resistance,  expression  of 
ERa  is  retained,  suggesting  that  resistance  involves 
either  changed  functionality  or  bypass  of  the 
receptor.  Culturing  breast  cancer  cells  in  estrogen- 
low  conditions  to  produce  long-term  estrogen 
deprivation  (LTED)  has  identified  mechanisms  of 
estrogen  hypersensitivity  and  estrogen  supersensitiv¬ 
ity  (Yue  et  al.  2002,  Martin  et  al.  2003,  2005a, b, 
Santen  et  al.  2005).  Estrogen  hypersensitivity  is 
characterized  by  the  ability  of  cells  to  respond  to 
levels  of  estrogen  at  concentrations  2-3  log  lower 
than  required  to  stimulate  wild- type  cells  (Yue  et  al. 

2002,  Santen  et  al.  2005).  This  mechanism  involves 
increased  expression  of  ERa  alongside  enhanced 
phosphorylation  of  ERa  Ser118  and  is  associated 
with  activation  of  the  ERK1/2  and  PI3-K  pathways. 
Estrogen  supersensitivity,  wherein  cells  are  appa¬ 
rently  estrogen  independent,  is  a  mechanism  again 
associated  with  enhanced  ERa  expression,  ERK 
activation,  and  activation  of  ERa  Ser118  and  involves 
ERa  being  supersensitized  by  growth  factor  acti¬ 
vation  (Martin  et  al.  2003,  2005a). 


While  higher  levels  of  ERa  expression  are 
generally  associated  with  enhanced  estrogen 
response,  in  certain  cases  tumors  expressing  high 
levels  of  ERa  can  be  insensitive  to  endocrine 
manipulation.  High  levels  of  ERa  expression  have 
been  associated  with  increased  proliferation  rates 
(Black  et  al.  1983)  and  poor  prognosis  in  breast 
cancer  patients  not  receiving  adjuvant  therapy 
(Black  et  al.  1983,  Thorpe  et  al.  1993).  It  has 
been  suggested  that  a  high  level  of  ERa  may  lead 
to  constitutive  activation  (Fowler  et  al.  2004).  This 
mechanism  has  recently  been  demonstrated  by 
Fowler  et  al.  (2004,  2006)  in  a  tetracycline- 
inducible  ERa  expression  model  of  the  MCF-7 
cell  line,  wherein  increased  ERa  expression 
resulted  in  aberrant  promoter  occupancy  and  gene 
activation  in  the  absence  of  estrogen.  The  increased 
receptor  activity  required  the  amino-terminal 
domain  and  was  not  inhibited  by  tamoxifen, 
supporting  the  notion  of  AF-1  activation,  yet  was 
independent  of  Ser104/106  and  Ser118  phosphoryl¬ 
ation  (Fowler  et  al.  2004). 

In  these  models,  the  expression  of  ERa  is  still 
critical  to  the  response  and  it  has  been  suggested  that 
use  of  a  SERD  such  as  fulvestrant  (faslodex,  ICI 
182  780)  would  be  a  beneficial  strategy  once  resistance 
to  aromatase  inhibitors  has  developed  (Johnston  et  al. 
2005,  Martin  et  al.  2005 b).  A  number  of  laboratories 
are  developing  models  of  resistance  to  this  agent  to 
identify  strategies  that  might  be  tried  at  the  onset  of 
resistance  (Dowsett  et  al.  2005,  Howell  2005,  Johnston 
et  al.  2005,  Martin  et  al.  2005 b,  Nicholson  et  al.  2005, 
Normanno  et  al.  2005). 

We  have  investigated  two  MCF-7  cell  lines 
(MCF-7/LCC1  and  MCF-7/LCC9),  which  have 
acquired  estrogen  insensitivity  and  with  variable 
sensitivity  to  tamoxifen  and  fulvestrant  to  identify 
novel  mechanisms  of  endocrine  resistance  that  might 
arise  in  clinical  specimens.  The  wild-type  ER-positive 
MCF-7  breast  cancer  cell  line  is  both  estrogen 
dependent  and  responsive  to  anti-estrogens,  such  as 
tamoxifen  and  fulvestrant.  The  MCF-7/LCC1  (LCC1) 
cell  line  was  derived  from  an  MCF-7  xenograft,  which 
had  grown  in  a  low  estrogen  environment  in  an 
immuno-deprived  mouse  and  which  was  known  to  be 
estrogen  independent  but  with  a  degree  of  estrogen 
sensitivity  (Brunner  et  al.  1993).  Treatment  of  the  cell 
line  with  fulvestrant  produced  the  MCF-7/LCC9 
(LCC9)  cell  line  which  is  fully  resistant  to  both 
estrogen  and  fulvestrant  (Brunner  et  al.  1997).  A 
number  of  novel  features  of  these  lines  were  identified 
within  this  study  and  are  reported  here. 
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Materials  and  methods 

Cell  proliferation 

MCF-7  cells  were  routinely  grown  in  phenol  red 
containing  Dulbecco’s  modified  Eagle  medium 
(DMEM)  supplemented  with  10%  fetal  calf  serum 
(FCS),  penicillin  (100  units/ml),  and  streptomycin 
(100  (g/ml).  LCC1  and  LCC9  cells  (source:  Dr  Robert 
Clarke,  V  T  Lombardi  Cancer  Research  Center, 
Georgetown  University  Medical  School,  Washington, 
DC,  USA)  were  routinely  kept  in  phenol-free 
containing  DMEM  supplemented  with  5%  dextran- 
activated  charcoal- stripped  fetal  calf  serum  (DCC), 
penicillin  (100  units/ml),  streptomycin  (100  (g/ml), 
and  2  mM  glutamine.  All  cells  were  grown  at  37  °C 
in  5%  C02.  To  determine  the  effects  of  17(3-estradiol 
(E2)  and  tamoxifen  on  cell  proliferation,  MCF-7  cells 
were  seeded  in  six-well  plates  in  phenol  red  containing 
DMEM  with  10%  fetal  bovine  serum  (FBS)  for  24  h. 
The  media  were  changed  to  phenol  red-free  DMEM 
with  5%  DCC  for  48  h.  The  cells  were  then 
supplemented  with  media  containing  either  1  nM  E2, 
1  pM  tamoxifen  or  both.  LCC1  and  LCC9  cells  were 
seeded  in  six- well  plates  in  phenol  red-free  containing 
DMEM  with  5%  DCC  and  after  24  h  supplemented 
with  E2  and/or  tamoxifen.  Cell  growth  was  evaluated 
using  a  Coulter  counter.  Fulvestrant  was  a  kind  gift 
from  Dr  Alan  Wakeling  (AstraZeneca,  Macclesfield, 
Cheshire,  UK).  For  studies  exploring  growth  in 
DMEM  without  serum,  the  sulforhodamine-B  (SRB) 
colorimetric  assay  was  used. 

Briefly,  log  phase  cells  were  seeded  into  96-well 
flatbottom  tissue  culture  plates.  The  following  day, 
cells  were  washed  in  PBS  and  media  replaced  with 
phenol  red-free  DMEM  for  48  h.  Cells  were  then 
treated  with  concentrations  of  E2  varying  from  10  fM 
to  1  pM  in  the  absence  or  presence  of  100  nM 
fulvestrant.  After  72  h,  plates  were  removed  from  the 
incubator  and  ice-cold  25%  trichloroacetic  acid  (TCA) 
solution  (50  pi)  added  to  each  well.  All  plates  were 
placed  on  ice  for  60  min  after  which  the  TCA  solution 
was  removed.  The  plates  were  washed  under  running 
water  and  dried  prior  to  staining  with  SRB  dye  solution 
(30  min  at  room  temperature)  and  the  trays  were 
washed  with  1%  glacial  acetic  acid  (X4)  at  room 
temperature,  air-dried,  and  resuspended  in  10  mM  Tris 
buffer  (pH  10.5;  150  pi)  before  reading  at  540  nm. 

RNA  extraction  and  RT-PCR 

Extraction  of  total  RNA  from  whole  cells  was 
performed  using  Tri-Reagent  (Sigma)  as  per  the 
manufacturers’  instructions.  RNA  concentration  was 


measured  using  a  spectrophotometer.  QuantiTect 
SYBR  Green  system  (Qiagen,  cat  no.  204243)  was 
used  according  to  the  manufacturers  instructions  for 
one  step  RT-PCR  in  a  total  of  15  pi  reaction  volumes, 
including  0.5  pM  each  primer  and  40  ng  RNA.  Real¬ 
time  cycler  conditions  were  RT:  50  °C  for  30  min; 
PCR:  initial  activation  95  °C  for  15  min  followed  by  40 
cycles  of  denaturation  94  °C  for  15  s,  annealing  57  °C 
for  30  s,  extension  72  °C  for  30  s,  and  a  final  extension 
of  72  °C  for  60  s.  The  following  primers  were  used: 

TFF1:  fwd  TTGTGGTTTTCCTGGTGTCA 
rev  CCGAGCTCTGGGACTAATCA 
ERa:  fwd  CCACCAACCAGTGCACCATT 
rev  GTCTTTCCGTATCCCACCTTTC 
PGR:  fwd  GTCAGTGGGCAGATGCTGTA 
rev  AGCCCTTCCAAAGGAATT 
ACTIN:  fwd  CTACGTCGCCCTGGACTTCGAGC 
rev  GATGGAGCCGCCGATCCACACGG 


Western  analysis 

Cells  were  washed  twice  with  PBS  and  lysed  in 
ice-cold  lysis  buffer  (50  mM  Tris  (pH  7.5),  5  mM 
EDTA  (pH  8.5),  150  mM  NaCl,  1%  Triton  X-100, 
aprotinin  10  pg/ml,  and  1  X  protease  cocktail  inhibitor 
(Roche)  for  10  min  and  the  debris  was  cleared  by 
centrifugation  at  13  000r.p.m.  for  6  min  at  4  °C). 
Protein  lysates  (100  pg)  were  resolved  on  7.5-12% 
SDS-PAGE  and  electrophoretically  transferred  to 
Immobilon-P  membranes.  After  transfer,  membranes 
were  blocked  and  probed  with  primary  antibody 
overnight  at  4  °C.  Immunoreactive  bands  were 
detected  using  chemiluminescent  reagents  (ECL  or 
SuperLuminol)  and  photographic  paper  (Hyperfilm, 
Amersham).  The  following  antibodies  were  used:  ERa 
(F-10;  Santa  Cruz  Biotech,  Santa  Cruz,  CA,  USA  sc- 
8002),  PGR  (ab-8;  Neomarkers,  Stratech  Scientific 
Ltd,  Newmarket,  Suffolk,  UK  (MS-298)),  P-ERK1/2 
(1:1000,  Cell  Signaling,  New  England  Biolabs, 
Hitchin,  Herts,  UK  #9101),  phospho-Ser118  ERa 
(1:500,  Cell  Signaling  #2511),  phospho-Ser167  ERa 
(1:500,  Cell  Signaling  #2514),  and  actin  (1:120  000, 
CP01,  Calbiochem,  La  Jolla,  CA,  USA).  Integrated 
optical  density  absorbance  values  were  obtained  by 
densitometric  analysis  using  a  gel  scanner  and 
analyzed  by  ‘Lab works’  gel  analysis  software  (UVP 
Life  Sciences,  Cambridge,  UK). 

Chromatin  immunoprecipitation  assays  (ChIP) 

Cells  were  grown  to  85-90%  confluence  in  phenol  red- 
free  DMEM  with  5%  DCC  for  at  least  48  h.  Cells  were 
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cross-linked  with  1%  formaldehyde  (37  °C  for  10  min) 
at  10-min  interval  over  a  90-min  time  course. 
Unreacted  formaldehyde  was  quenched  by  gentle 
agitation  at  room  temperature  for  10  min  with 
0.125  M  glycine.  Cells  were  then  washed  twice  with 
ice-cold  PBS,  collected  into  PBS  containing  protease 
inhibitors  (Roche),  and  centrifuged  for  4  min  at 
2000  r.p.m.  at  4  °C.  The  pellets  were  resuspended  in 
lysis  buffer  (1%  SDS,  10  mM  EDTA,  50  mM  Tris-HCl 
(pH  8.1),  and  1  X  protease  inhibitor  cocktail), 
incubated  on  ice  for  10  min,  and  sonicated  (12X20  s 
at  two  amplitude  microns,  Soniprep  150,  MSE)  to 
fragment  DNA  to  ~500  bp.  Following  centrifugation 
for  15  min  at  13  000  r.p.m.  and  4  °C,  supernatants  were 
collected  and  resuspended  in  dilution  buffer  (0.01% 
SDS,  1%  Triton  X-100,  1.2  mM  EDTA,  16.7  mM 
Tris-HCl  (pH  8.1),  167  mM  NaCl,  and  IX  protease 
inhibitor  cocktail).  Chromatin  were  precleared  with 
1  pg  anti-rabbit  or  anti-mouse  IgG,  2  pg  sheared 
salmon  sperm  DNA,  and  Protein-G- Agarose  (50  pi  of 
50%  slurry  in  dilution  buffer)  for  3  h  at  4  °C. 
Immunoprecipitation  using  Protein-G-Agarose  Beads 
(Roche)  was  performed  overnight  at  4  °C  with  anti- 
ERa  HC-20  antibody  (sc-543,  Santa  Cruz).  Beads  were 
washed  sequentially  for  5  min  each  at  4  °C  with  TSE  I 
(20  mM  Tris  (pH  8.1),  2  mM  EDTA,  150  mM  NaCl, 
1%  Triton  X-100,  and  0.1%  SDS),  TSE  II  (20  mM  Tris 
(pH  8.1),  2  mM  EDTA,  500  mM  NaCl,  1%  Triton 
X-100,  and  0.1%  SDS),  and  buffer  III  (10  mM  Tris  (pH 
8.1),  0.25  M  LiCl,  1  mM  EDTA,  1%  NP40,  and  1% 
deoxycholate).  Precipitates  were  then  washed  twice 
with  TE  buffer  and  the  protein/DNA  complexes  were 
eluted  twice  with  0.1  M  NaHC03  and  1%  SDS.  Heat 
treatment  at  65  °C  overnight  reversed  formaldehyde 
cross-links.  DNA  fragments  were  purified  using 
QIAquick  Spin  Kit  columns  (Qiagen)  and  amplified 
using  the  QuantiTect  SYBR  Green  system  (Qiagen,  cat 
no.  204242).  TFF1  PCR  conditions  were:  initial 
activation  of  95  °C  for  15  min  followed  by  45  cycles 
of  94  °C  for  15  s,  55  °C  for  30  s,  72  °C  for  30  s,  and  a 
final  extension  of  72  °C  for  5  min.  TFF1  primer 
sequences:  fwd  GACGGAATGGGCTTCATGAGC 
and  rev  CTGAGACAATAATCTCCACTG.  For  the 
distal  region,  primers  were:  fwd  GAGTTTGGCCTCC- 
CACATTA  and  rev  CTTGCCTCTGCATTCTCTCC. 

Short  interfering  (siRNA)  transfections 

MCF-7  cells  were  seeded  at  0.5  X 106  cells  per  T75  flask 
in  DMEM  as  mentioned  previously.  After  24  h,  the 
media  were  changed  to  phenol  red-free  containing 
DMEM  with  5%  DCC  for  48  h.  LCC1  and  LCC9  cells 
were  seeded  directly  into  phenol  red-free  containing 


DMEM  with  5%  DCC  for  24  h  prior  to  transfection. 
Cells  were  transfected  with  siRNA  for  4  h  using 
Oligofectamine  reagent  (Invitrogen)  after  which  time 
1  nM  E2  was  added  for  a  further  48  h  prior  to  RNA  and 
protein  extraction.  For  the  7-day  time  course,  the  media 
were  left  unchanged  after  the  initial  changes.  For  siRNA 
growth  assays,  cells  were  seeded  as  for  growth 
characterization  as  mentioned  previously.  siRNA 
transfections  were  carried  out  as  described  earlier  but 
scaled  down  for  24- well  plates.  Following  siRNA 
treatment  for  4  h,  cells  were  treated  with  1  nM  E2  or 
100  nM  fulvestrant  or  a  combination  and  cell  counts  on 
days  0,  3,  and  6  were  estimated  using  a  Coulter  counter. 
The  following  siRNA  sequences  were  used:  ER  RNAi 
1;  ESR1  SMARTpool  (four  pooled  sequences;  Upstate 
Biotechnology,  Lake  Placed,  NY,  USA;  M-003401; 
100  nmol),  ER  RNAi  2;  5 A  A  AC  AGG  AGG  A  A- 
GAGCTGCCA  (Ambion;  40  nmol),  ER  RNAi  3;  5'- 
AACCTCGGGCTGTGCTCTTTT  (Ambion,  Hunting¬ 
don,  Cambridgeshire,  UK;  40  nmol),  and  negative 
RNAi:  Upstate  (D-001206;  100  nmol). 

Results 

Increased  ERa  expression  in  resistant  cell  lines 

To  explore  the  possibility  that  high  ERa  expression 
leads  to  estrogen-independent  growth  in  endocrine- 
resistant  cells,  the  expression  levels  of  ERa  in  resistant 
lines  (LCC1  and  LCC9)  were  compared  with  levels  in 
wild-type  MCF-7  cells.  Both  resistant  lines  expressed 
between  four-  and  elevenfold  more  ERa  mRNA  than 
wild- type  cells  (Fig.  1A).  ERa  protein  levels  were 
clearly  elevated  in  LCC1  cells  relative  to  MCF-7  cells 
(sevenfold)  and  less  markedly  in  LCC9  cells  (Fig.  IB). 
E2  decreased  ERa  protein  in  MCF-7  cells  at  48  h  and 
this  has  been  explained  by  proteosomal  degradation,  a 
process  speculated  to  limit  the  action  of  estrogen 
signaling  (Nawaz  et  al.  1999;  Fig.  1C).  Similarly,  both 
resistant  lines  demonstrated  ERa  turnover,  suggesting 
that  ERa  is  binding  to  E2  in  all  cases.  In  contrast, 
tamoxifen  treatment  results  in  maintenance  of  the 
receptor  expression  levels  in  all  three  cell  lines 
(Fig.  1C). 

Addition  of  1  nM  17  (3-estradiol  (E2)  to  MCF-7  cells 
produced  a  marked  stimulation  of  growth  to  cells 
cultured  in  estrogen-depleted  (double  charcoal- 
stripped  FCS)  medium  (Fig.  2A).  In  the  absence  of 
E2,  MCF-7  cells  are  essentially  static  (Fig.  2A).  In 
contrast,  LCC1  cells  grow  rapidly  in  estrogen-depleted 
conditions  and  show  an  approximately  twofold 
stimulation  of  growth  on  addition  of  E2  (Fig.  2B). 
LCC9  cells  showed  a  lack  of  response  to  E2,  again 
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Figure  1  ERa  expression  in  MCF-7,  LCC1,  and  LCC9  cells.  (A)  ERa  mRNA  expression.  Cells  were  grown  in  charcoal-stripped 
serum-containing  medium  for  at  least  48  h  and  RNA  was  collected.  A  representative  experiment  is  shown  of  at  least  two  experiments 
carried  out.  Each  column  presents  mean  of  triplicate  RT-PCR  analysis  for  each  sample  demonstrating  mRNA  expression  relative  to 
actin  expression.  Error  bars  =  s.D.  Statistical  significance  noted  for  treatment  groups  versus  matched  control  (one-way  ANOVA  and 
multiple  comparison  Tukey-Kramer  test;  *P<0.05).  (B)  Western  blot  analysis  of  ERa  (66  kDa)  in  breast  cancer  cell  lines  grown  in 
charcoal-stripped  serum-containing  medium  for  48  h  prior  to  protein  collection.  One  hundred  micrograms  of  protein  were  loaded  per 
lane  and  detected  using  either  anti-ERa  (Santa  Cruz  Biotech)  or  anti-actin  (Calbiochem)  antibodies  as  described  in  Materials  and 
methods.  (C)  Western  blot  analysis  of  ERa  (66  kDa)  in  breast  cancer  cell  lines  grown  in  charcoal-stripped  serum-containing  medium 
for  at  least  48  h  prior  to  protein  collection.  One  hundred  micrograms  of  protein  were  loaded  per  lane  and  detected  using  either  anti- 
ERa  (Santa  Cruz  Biotech)  or  anti-actin  (Calbiochem)  antibodies  as  described  in  Materials  and  methods.  (D)  ERa  binding  to  the  TFF1 
promoter.  Basal  recruitment  of  ERa  to  the  TFF1  promoter  was  determined  by  ChIP  analysis  on  untreated  cells.  The  ChIP  method 
used  was  as  described  in  Materials  and  methods  and  immunoprecipitated  TFF1  promoter  was  quantified  by  real-time  PCR.  The 
input-corrected  expression  values  were  determined  by  normalizing  to  the  inputs.  Data  are  presented  as  mean  +  s.E.  Groups  were 
compared  with  the  Kruskal-Wallis  test  with  Dunn’s  multiple  comparison  test  (*P<0.05).  Binding  to  the  promoter  region  is  compared 
with  binding  to  a  region  3.5  kb  distal  to  the  promoter  wherein  only  background  binding  was  observed. 


growing  very  rapidly  in  the  absence  of  E2  (Fig.  2C). 
Addition  of  1  pM  tamoxifen  to  MCF-7  cells  antag¬ 
onized  the  E2-stimulated  growth  in  this  cell  line. 
Tamoxifen  also  inhibited  the  E2-stimulated  growth  of 
LCC1  cells  but  had  no  effect  on  LCC9  cells  (Fig.  2B 
and  C).  These  results  are  consistent  with  wild- type 
cells  being  estrogen  dependent,  LCC1  cells  demon¬ 
strating  partial  estrogen  dependence  and  LCC9  cells 
being  fully  estrogen  independent. 

Reduced  ERa  Ser118  phosphorylation 
in  LCC9  cells 

Several  frequently  cited  mechanisms  of  estrogen- 
independent  activation  of  ERa  involve  phosphoryl¬ 
ation  of  ERa  at  the  Ser118  or  Ser167  residues  mediated 
via  ERK  or  Akt  respectively  (Bunone  et  al.  1996, 
Martin  et  al.  2000,  Lannigan  2003).  While  the  Ser118 


residue  is  a  major  site  of  E2-induced  phosphorylation, 
Ser167  is  not  (Lannigan  2003).  The  latter  site  is 
activated  by  growth  factor  signaling.  In  view  of  these 
previous  observations,  we  first  investigated  whether 
ERa  Ser118  or  Ser167  phosphorylation  were  increased 
in  the  absence  of  estrogen  in  the  resistant  cell  lines. 
Neither  was  there  evidence  of  increased  Ser118 
phosphorylation  in  the  resistant  lines  relative  to 
MCF-7  under  basal  conditions,  nor  was  Ser167 
phosphorylation  increased  (Fig.  3A-C).  Furthermore, 
phospho-ERKl/2  expression  was  unchanged  in  the 
lines  (Fig.  3C).  On  E2  addition,  there  was  a  marked 
increase  in  Ser118  phosphorylation  in  MCF-7  cells  and 
this  was  also  observed  in  the  LCC1  cell  line  (Fig.  3 A 
and  B).  However,  minimal  change  was  observed  on  E2 
addition  to  LCC9  cells  (Fig.  3 A  and  B).  Ser118 
phosphorylation  has  been  proposed  to  affect  cofactor 
recruitment  and  this  might  explain  the  reduced 
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Figure  2  Growth  characterization  of  MCF-7  and  MCF-7  variant  cells.  (A)  MCF-7  cells,  (B)  LCC1 ,  and  (C)  LCC9  cells  were  plated  for 
24  h  and  maintained  in  reduced  media  for  48  h  before  treatment.  Cells  were  then  left  untreated  (control  group),  treated  with  1  nM  E2, 
1  |llM  tamoxifen  or  1  nM  E2  and  1  tamoxifen.  Cells  were  counted  on  day  0  (72  h  after  plating/day  of  treatment  start)  and  days  2/4/6 
using  a  Coulter  counter.  Mean  cell  counts  of  triplicate  samples  and  duplicate  counts  for  each  time  point  in  each  treatment  group  are 
expressed.  Error  bars  =  s.D.  A  representative  experiment  is  shown  of  at  least  four  experiments  carried  out. 


transcriptional  (as  mentioned  below)  and  growth 
responses  observed  on  E2  addition  to  this  cell  line. 
Tamoxifen  alone  produced  a  small  increase  in  Ser118 
phosphorylation  in  MCF-7  and  LCC1  cells  but  not  in 
LCC9  cells  (Fig.  3 A  and  B).  Tamoxifen  also  produced 
a  reduction  of  estrogen’s  Ser118  phosphorylation  in  the 
MCF-7  and  LCC1  cell  lines  (Fig.  3 A  and  B). 

Modified  DNA  binding  of  ERa  in  resistant 
cell  lines 

To  explore  whether  high  ERa  expression  was  reflected  in 
enhanced  DNA  binding  in  the  absence  of  E2,  ChIP 
methodology  was  used  to  examine  ERa  binding  to  the 
promoter  of  the  E2-responsive  gene  TFF1  in  the  MCF-7, 
LCC1,  and  LCC9  cell  lines.  LCC9  cells  had  >  2.5-fold 
greater  ERa  binding  to  the  TFF1  promoter  than  MCF-7 
cells  (Fig.  ID).  However,  this  binding  was  significantly 
higher  in  LCC1  cells  with  levels  greater  than  eightfold 
above  MCF-7  cells.  This  enhanced  ERa  binding  in  LCC1 
cells  was  equivalent  to  the  increased  expression  of  ERa 
protein  and  is  consistent  with  the  suggestion  by  Fowler  et 
al.  (2004)  that  enhanced  ERa  protein  expression  can  lead 
to  increased  DNA  binding.  As  a  control,  binding  to  a 
region  3.5  kb  distal  to  this  region  indicated  only 
background  levels  as  expected  (Fig.  ID). 


Growth  responses  to  estrogen  and  tamoxifen 
in  the  wild-type  and  variant  cell  lines  are  reflected 
in  transcriptional  changes 

To  investigate  the  differences  in  estrogen  and  anti¬ 
estrogen  activation  processes,  indicator  genes  that 
reflected  the  different  growth  responses  were  next 
investigated.  Transcriptional  changes  in  the  estrogen- 
regulated  genes  TFF1  and  PGR  were  measured 
and  modulated  expression  was  compared  with  the 
growth  changes. 

Expression  of  TFF1  mRNA  in  the  absence  of  E2  was 
higher  in  both  resistant  lines  compared  with  MCF-7  cells 
(Fig.  4A).  After  48-h  E2  (1  nM)  treatment,  TFF1  mRNA 
was  increased  by  >  20- fold  in  MCF-7  cells,  but  only  one- 
to  twofold  in  the  resistant  lines  although  this  increase  was 
significant.  Tamoxifen  (1  jliM)  produced  a  small  increase 
in  TFF1  expression  in  MCF-7  and  LCC1  cells  but  not  in 
the  LCC9  cell  line  (Fig.  4A).  These  levels  broadly  reflect 
the  growth  differences  observed. 

The  expression  of  PGR  mRNA  in  the  absence  of  E2 
was  greater  in  LCC1  and  LCC9  lines  compared  with 
MCF-7  cells  (Fig.  4B).  As  for  TFF1,  after  48-h  E2 
treatment,  PGR  mRNA  was  increased  by  >  20-fold  in 
MCF-7  cells  and  2-5-fold  in  LCC1  and  LCC9  cell  lines 
(Fig.  4B).  Tamoxifen  also  increased  the  PGR  mRNA 
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Figure  3  Effects  of  E2  and  tamoxifen  on  ERa  phosphorylation  in  the  resistant  cell  lines.  (A)  Western  analysis  of  ERa  phospho-Ser118 
after  30-min  treatment  with  control  (no  treatment),  1  nM  E2,  1  riM  tamoxifen,  or  1  nM  E2  and  1  jiM  tamoxifen.  Lysates  were  run  on  a 
10%  SDS  gel  and  membranes  probed  with  anti-phospho  ERa  Ser118  antibody  (1:1000).  Lysates  were  also  probed  for  actin 
expression  to  compare  protein  loading.  (B)  Histogram  representing  optical  densities  from  triplicate  western  blots  of  ERa  phospho- 
Ser1 18  after  30-min  treatment  with  control  (no  treatment),  1  nM  E2, 1  jiM  tamoxifen,  or  1  nM  E2and  1  |iM  tamoxifen.  Values  were  actin 
corrected  and  then  standardized  to  a  control  sample.  The  control  sample  was  a  30  min  Entreated  MCF-7  sample  and  was  used  on 
all  gels  as  a  standard  to  allow  comparisons  between  runs.  Statistical  comparisons  were  made  with  each  cell  line’s  control  level; 
*P<0.01;  fP<  0.001  (ANOVA).  (C)  Western  analysis  of  ERa  phospho-Ser167  and  phospho-ERK1/2  in  cell  lines.  Untreated  lysates 
were  probed  with  antibodies  specific  for  ERa  phospho-Ser167  and  phospho-ERK1/2.  A  positive  control  lane  of  MCF-7  cells  treated 
with  1  nM  HRGp  was  used.  Lysates  were  also  probed  for  actin  expression. 


expression  level  not  only  in  MCF-7  cells,  but  also  in 
LCC1  cells  producing  effects  equivalent  to  that  of  E2  in 
the  latter  cell  line.  No  change  was  observed  in  the 
LCC9  cell  line. 

These  results  are  consistent  with  transcription  of 
TFF1  and  PGR  being  increased  by  ligand-independent 
mechanisms  in  LCC1  and  LCC9  cell  lines  with 
estrogen  and  tamoxifen  producing  an  additional 
ligand-dependent  increase. 

Effect  of  removal  of  ERa  on  the  growth 
of  the  cell  lines 

To  determine  the  relative  importance  of  ERa  on 
downstream  gene  expression  and  growth  of  MCF-7, 
LCC1,  and  LCC9  cells,  we  investigated  the  effects  of 
removing  ERa,  either  by  specific  siRNA  inhibition  of 
receptor  synthesis  or  through  inhibition  and 
degradation  of  the  receptor  by  fulvestrant. 

A  panel  of  interfering  RNAs  (siRNAs)  were  initially 
compared  for  their  ability  to  transiently  reduce  ERa 
expression  and  were  transfected  into  the  MCF-7  cell 


line.  RNAi  1  is  a  pooled  set  of  four  targeted  sequences 
(Imai  et  al.  2005)  while  RNAi  2  ^'-AAACAGGAG- 
GAAGAGCTGCCA)  and  RNAi  3  (S'-AACCT- 
CGGGCTGTGCTCTTTT)  are  individually  targeted 
sequences  (Leu  et  al  2004).  Of  the  three,  RNAi  2 
produced  the  best  reduction  of  ERa  mRNA  and  protein 
and  was  selected  for  further  experiments  (Fig.  5A  and 
B).  Quantitative  RT-PCR  analysis  showed  that,  48  h 
after  transfection,  ERa  RNAi  2  treatment  resulted  in  an 
85%  decrease  in  ERa  mRNA  expression  and  an  87% 
decrease  in  the  presence  of  E2  (Fig.  5C).  LCC1  and 
LCC9  cells  have  significantly  higher  basal  expression 
of  ERa  mRNA  and  siRNA  removal  caused  an  82  and 
73%  decrease  respectively  with  similar  reductions  in 
the  presence  of  E2  (Fig.  5C).  Western  analysis  of  the 
MCF-7  and  LCC1  cell  lines  demonstrated  that  RNAi  2 
produced  ERa  protein  knockdown  over  a  7-day  period 
(Fig.  5D)  and  it  was  effective  in  all  three  cell  lines 
(Fig.  5E).  This  reduction  in  ERa  protein  was 
accompanied  by  a  decrease  in  PGR  protein  (Fig.  5E). 
Thus,  it  appeared  that  gene  expression  in  all  three  cell 
lines  was  ERa  dependent. 
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Figure  4  Effect  of  estrogen,  tamoxifen,  and  fulvestrant  on  (A)  TFF1  expression  and  (B)  PGR  mRNA  expression  in  the  cell 
lines.  Relative  mRNA  expression  values  of  TFF1  and  PGR  in  the  cell  lines  were  measured  by  real-time  RT-PCR  using 
specific  primer  pairs.  RNA  was  collected  at  48  h  and  was  extracted  from  either  untreated  (control)  cells  or  cells  treated  with 
1  nM  E2,  1  itM  tamoxifen,  1  nM  E2  +  1  itM  tamoxifen,  100  nM  fulvestrant,  or  1  nM  E2  +  100nM  fulvestrant.  Each  column 
represents  mean  of  triplicate  PCR  analysis  for  each  sample  demonstrating  mRNA  expression  relative  to  actin  expression. 
Error  bars  =  s.D.  Statistical  significance  noted  for  treatment  groups  compared  to  matched  control  where  *P< 0.05,  untreated 
control  versus  treatment  group;  tP<0.05,  E2  control  versus  treatment  group  (ANOVA  and  multiple  Tukey-Kramer 
comparison  test). 


This  was  investigated  further  using  fulvestrant. 
Fulvestrant  abrogates  E2-induced  gene  transcription 
by  binding,  blocking,  and  causing  the  degradation  of 
ERa  (Parker  1993).  Fulvestrant  treatment  in  MCF-7 
cells  blocked  E2-induced  expression  of  TFF1  and  PGR 
(Fig.  4 A  and  B).  In  addition,  ligand-independent  and 
E2-induced  TFF1  and  PGR  expression  in  LCC1  cells 
were  reduced  on  fulvestrant  treatment.  These  data 
confirm  that  for  LCC1  cells  TFF1  and  PGR  induction 
are  dependent  on  ERa  expression.  However,  LCC9 
cells  are  resistant  to  fulvestrant  treatment  and  as  such 
no  change  in  TFF1  expression  and  only  a  minor  change 
in  PGR  expression  was  observed.  The  effect  of 
fulvestrant  on  the  growth  of  all  three  cell  lines  was 
also  investigated  in  the  complete  absence  of  serum 
(Fig.  6).  Under  these  conditions,  MCF-7  cells  did  not 
grow  over  a  72-h  period.  LCC1  cells,  however,  still 
proliferated  and  the  addition  of  E2  had  little  effect  on 
growth  confirming  their  independence  of  E2.  Under 
these  conditions,  fulvestrant  was  able  to  oppose  the 
effect  of  low  concentrations  of  E2  again  indicating 
dependence  on  ERa.  In  contrast,  LCC9  cells 
were  completely  insensitive  to  both  E2  and 


fulvestrant.  Fulvestrant  degraded  ERa  protein  in  all 
three  lines  which  is  shown  in  Fig.  7 A. 

To  determine  how  critical  levels  of  ERa  expression 
were  for  the  growth  of  MCF-7,  LCC1,  and  LCC9  cell 
lines,  we  used  RNAi  removal  with  or  without  fulvestrant 
to  inhibit  the  synthesis  of  ERa  protein  (Fig.  7B-D).  E2- 
induced  MCF-7  cell  growth  was  significantly  decreased 
(33%)  by  ERa  removal  and  abolished  by  all  com¬ 
binations  of  fulvestrant  alone  or  with  RNAi.  LCC1  cells 
grew  in  the  absence  of  E2  and  RNAi  removal  had  only  a 
minor  effect  on  growth.  E2-induced  LCC1  cell  growth 
was  reduced  by  approximately  40%  when  ERa  was 
removed  through  RNAi,  but,  unlike  MCF-7  cells, 
fulvestrant  alone  was  not  enough  to  abolish  growth  - 
this,  however,  could  be  accomplished  though  through 
combination  with  RNAi.  LCC9  cell  growth  in  the 
absence  of  E2  was  reduced  by  ERa  RNAi.  A  similar 
decrease  was  observed  in  the  presence  of  E2.  LCC9  cells 
are  fulvestrant  resistant  and  no  effect  on  growth  was 
observed  with  this  agent.  No  combinations  of  fulvestrant 
or  RNAi  were  able  to  totally  abolish  growth.  These 
results  indicate  a  varying  degree  of  dependence  on  ERa 
for  growth  in  the  three  cell  lines. 
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Figure  5  Effects  of  ERa  RNAi  on  ERa  and  PGR  expression  in  the  cell  lines.  (A)  Expression  of  ERa  mRNA  after  treatment  with  a 
range  of  ERa  mRNA-targeted  RNAis.  ERa  mRNA  expression  was  measured  by  quantitative  RT-PCR  of  mRNA  from  MCF-7  cells 
48  h  after  RNAi  treatment  in  the  presence  of  1  nM  E2.  Data  are  presented  as  mean  +  s.D.  of  actin-corrected  values  from  triplicate 
samples.  The  RNAi  transfection  method  and  RNAi  sequences  used  are  described  in  Materials  and  methods.  Statistical  significance 
noted  for  treatment  groups  compared  with  matched  control  where  *P<  0.05,  untreated  control  versus  treatment  group  (ANOVA  and 
multiple  Tukey-Kramer  comparison  test).  (B)  Western  analysis  of  ERa  protein  expression  in  MCF-7  cells  48  h  after  siRNA  treatment. 
ERa  was  probed  with  the  F-10  antibody  and  actin  is  shown  as  a  loading  control.  (C)  Expression  of  ERa  mRNA  after  treatment  with 
RNAi  2.  ERa  mRNA  expression  was  measured  by  quantitative  RT-PCR  of  mRNA  from  cell  lines  48  h  after  RNAi  2  treatment  in  the 
presence  or  absence  of  1  nM  E2.  Data  are  presented  as  mean  +  s.D.  of  actin-corrected  values  from  triplicate  samples.  The  RNAi 
transfection  method  and  RNAi  sequences  used  are  described  in  Materials  and  Methods.  Statistically  significance  differences  are 
noted  for  treatment  groups  compared  with  matched  control  where  *P<  0.05,  untreated  control  versus  treatment  group  (ANOVA  and 
multiple  Tukey-Kramer  comparison  test).  (D)  Western  analysis  time  course  of  the  effect  of  RNAi  2  treatment  on  ERa  protein 
expression  in  the  MCF-7  and  LCC1  cell  lines.  ERa  was  probed  with  the  F-10  antibody  and  actin  is  shown  as  a  loading  control.  (E) 
Western  analysis  of  ERa  and  PGR  protein  expression  in  the  cell  lines  48  h  after  siRNA  treatment.  ERa  was  probed  with  the  F-10 
antibody,  PGR  with  Ab  8  and  actin  is  shown  as  a  loading  control. 


Discussion 

Aromatase  inhibitors  are  now  used  for  the  adjuvant 
treatment  of  most  hormone  receptor-positive  early 
breast  cancer.  Despite  the  improvement  they  offer  over 
tamoxifen  alone,  recurrences  still  occur,  and  thus 
models  of  resistance  to  both  tamoxifen  and  estrogen 
deprivation  are  required.  The  series  of  MCF-7 -derived 
cell  lines  provides  an  excellent  model  system  for  the 
exploration  of  mechanisms  of  stepwise  acquisition  of 
resistance  to  tamoxifen  and  estrogen  deprivation.  Most 
models  to  date  have  been  derived  in  vitro ,  which  makes 
LCC1  cells  interesting  as  the  initial  estrogen  depri¬ 
vation  was  achieved  in  vivo  and  therefore  might  reflect 
features  that  could  arise  in  a  primary  breast  cancer 
(Brunner  et  al.  1993).  In  many  of  the  in  vitro- derived 
LTED  models,  acquired  resistance  is  due  to  enhanced 
sensitization  to  low  concentrations  of  estrogen,  which 


often  involves  crosstalk  with  growth  factor- signaling 
pathways  (Martin  et  al.  2003,  2005a, b).  LCC1  cells 
have  certain  of  the  characteristics  of  the  LTED 
phenotype  (Yue  et  al.  2002,  Martin  et  al.  2003, 
2005a, b,  Santen  et  al.  2005)  such  as  a  higher 
expression  level  of  ERa,  an  ability  to  grow  in  low- 
estrogen  conditions  and  elevated  TFF1  expression.  The 
continuous  culturing  of  LCC1  cells  in  low  estrogen 
conditions  may  well  contribute  to  the  increased 
expression  of  ERa  in  this  cell  line. 

However,  unlike  most  LTED-derived  cells,  which 
show  little  response  to  physiological  levels  of  estrogen 
yet  are  sensitive  to  very  low  levels  of  estrogen,  LCC1 
cells  appear  truly  insensitive  to  the  addition  of  low 
levels  of  exogenous  estrogen.  Similarly,  while  most 
LTED  cells  show  basal  activation  of  ERK1/2  acti¬ 
vation  and  ERa  via  Ser118  phosphorylation,  LCC1  and 
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Figure  6  Effect  of  the  growth  of  the  cell  lines  in  serum-free 
media  and  treated  with  varying  concentrations  of  E2  in  the 
absence  or  presence  of  fulvestrant  (1  nM).  Cells  were  plated 
and  after  establishment  placed  in  serum-free  medium  for  48  h. 
E2  with  or  without  fulvestrant  was  added  and  plates  left  for  72  h. 
Relative  cell  numbers  were  then  assessed  by  SRB  assay  as 
described  in  Materials  and  methods.  Inset  in  MCF-7  figure: 
effect  of  E2  on  MCF-7  cells  grown  in  5%  double  charcoal- 
stripped  fetal  serum. 

LCC9  cells  do  not.  The  ER,  however,  is  still  clearly 
functional  in  LCC1  cells  and  linked  to  growth 
regulation  as  estrogen  addition  can  produce  an  increase 
in  growth  which  could  be  reversed  by  tamoxifen.  ERa 
is  also  downregulated  by  the  addition  of  estrogen  and 
markedly  phosphorylated  at  Ser118.  Additionally,  the 
ERa  downregulator  fulvestrant  reduces  expression  of 
TFF1  and  inhibits  growth.  These  effects  are  more 
marked  when  cells  are  exposed  to  fulvestrant  with 
siRNA  removal  of  ERa. 

While  constitutive  activation  of  ERa  may  be 
achieved  in  some  instances  by  phosphorylation  of 


Ser118  mediated  by  growth  factor-driven  activation  of 
ERK,  an  increased  expression  of  ERa  alone  might 
account  for  increased  DNA  binding.  In  support  of  this, 
there  was  enhanced  binding  of  ERa  to  the  TFF1 
promoter  in  the  absence  of  added  estrogen  in  both  the 
LCC1  and  LCC9  cell  lines.  In  addition,  TFF1 
transcription  was  markedly  increased  in  the  resistant 
cell  lines  consistent  with  this  enhanced  ERa-binding 
driving  transcription.  Direct  support  for  such  a 
mechanism  has  recently  been  demonstrated  in  an 
MCF-7  cell  line  using  a  tetracycline-inducible  ERa 
overexpression  model  (Fowler  et  al.  2004,  2006).  As 
with  the  data  mentioned  earlier,  the  results  suggested 
that  elevated  levels  of  ERa  resulted  in  activation  of 
receptor  transcriptional  function  in  a  manner  distinct 
from  mechanisms  that  involve  ligand  binding  or 
growth  factor-induced  phosphorylation  of  the  Ser104, 
Ser106  or  Ser118  sites.  The  mechanism  required  the 
amino-terminal  A/B  domain  and  was  not  inhibited  by 
tamoxifen.  It  was  also  uncoupled  from  ERK  activation. 
The  hypothesis  proposed  was  that  overexpression  of 
unliganded  ERa  stabilized  interactions  with  the  basal 
transcriptional  machinery,  which  at  normal  receptor 
levels  may  be  too  weak  to  support  effective  transcrip¬ 
tion  (Fowler  et  al.  2004). 

These  results  together  support  a  model  wherein 
growth  (and  TFF1  transcriptional  activation)  in  LCC1 
cells  is  dependent  on  ERa.  This  dependency  has  some 
ligand  (i.e.,  estrogen)  responsiveness  but  is  largely 
ligand  independent.  The  ligand-dependent  component 
may  be  reversed  by  tamoxifen.  The  ligand  indepen¬ 
dency  appears  to  involve  neither  growth  factor 
activation  via  the  Ser118  or  Ser167  phosphorylation 
routes  nor  hypersensitization  (where  low  levels  of 
estrogen  produce  apparent  independence).  Instead  the 
ligand  independence  appears  to  be  explained  by  the 
high  level  of  ER  expression  leading  to  constitutive 
activation  and  promoting  DNA  binding  and  transcrip¬ 
tional  activation. 

We  have  shown  that  ERa  is  functionally  active  in 
the  LCC1  model  and  since  this  has  also  been  shown  in 
models  demonstrating  LTED,  a  logical  clinical 
strategy  to  attempt  after  development  of  resistance  in 
a  low  estrogen  environment  (such  as  produced  by 
aromatase  inhibitor  treatment)  is  to  downregulate  the 
receptor  using  fulvestrant  (Johnston  et  al.  2005,  Martin 
et  al.  2005a, b).  This  strategy  clearly  is  effective  at 
inhibiting  growth  in  LCC1  cells.  However,  the  LCC9 
variant  was  derived  after  exposure  and  development  of 
resistance  to  fulvestrant  (Brunner  et  al.  1997)  and 
showed  no  growth  response  to  either  estrogen  or 
tamoxifen.  In  this  cell  line,  the  negligible  changes  of 


1130 


www.endocrinology-journals.org 


Endocrine-Related  Cancer  (2006)  1 3  1 1 21  -1 1 33 


B  MCF-7 


D 


Figure  7  Effect  of  fulvestrant  on  ERa  expression  and  combined  with  ERa  siRNA  on  the  growth  of  the  cell  lines.  (A)  Western  blot 
analysis  of  ERa  (66  kDa)  in  breast  cancer  cell  lines  in  control,  1  nM  E2,  100  nM  fulvestrant,  or  1  nM  E2  and  100  nM  fulvestrant- 
treated  groups  at  48  h.  One  hundred  micrograms  of  protein  were  loaded  per  lane  and  detected  using  anti-ERa  (Santa  Cruz  Biotech) 
antibody  as  described  in  Materials  and  Methods.  Actin  expression  is  also  shown.  (B-D)  Effects  of  fulvestrant,  ERa  siRNA,  or 
combinations  on  the  growth  of  the  cell  lines.  (B)  MCF-7,  (C)  LCC1 ,  and  (D)  LCC9  cells  were  treated  with  1  nM  E2, 1 00  nM  fulvestrant, 
40  nM  ERa  siRNA,  or  combinations  of  these.  siRNA  treatment  was  for  4  h  only,  while  E2  and  fulvestrant  were  present  throughout  the 
time  course.  Comparisons  are  made  with  the  negative  siRNA  control  which  gave  an  equivalent  growth  effect  to  no  treatment.  Data 
are  presented  as  mean  +  s.E.  from  quadruplicate  samples.  Statistical  significance  noted  for  treatment  groups  compared  with 
matched  control  where  *P<  0.05,  negative  RNAi  control  versus  treatment  group;  tP<0.05,  negative  RNAi  +  E2  control  versus 
treatment  group  (ANOVA  and  multiple  Tukey-Kramer  comparison  test). 


ERa  Ser118  phosphorylation  obtained  on  estrogen  or 
tamoxifen  addition  contrasted  with  observations  in  the 
other  cell  lines.  Markedly  reduced  phosphorylation  is 
likely  to  affect  cofactor  binding  and  our  initial  findings 
suggest  that  pi 60  binding  (specifically  AIB1)  is 
reduced  in  this  cell  line,  again  consistent  with 
endocrine  insensitivity  (Kuske  et  al.  2004).  However, 
it  is  quite  clear  that  fulvestrant  can  downregulate  the 
receptor  and  even  extremely  high  levels  of  fulvestrant 
(10  pM)  were  unable  to  influence  growth  (data  not 
shown).  Despite  this,  siRNA  removal  of  ERa  produced 
some  growth  inhibition  suggesting  a  reduced  but  still 
measurable  dependency  on  ERa. 

In  conclusion,  these  results  suggest  that  multiple 
changes  contribute  to  endocrine  resistance.  While  ER 
still  demonstrates  functionality  in  LCC1  cells,  there  is  a 
major  shift  to  ligand  independence.  This  independence 
can  be  explained  by  the  high  level  of  ER  expression 
found  in  these  cells  and  could  lead  to  constitutive 
activation  of  the  receptor.  These  cells  still  show  a  degree 
of  dependency  on  estrogen  and  this  can  be  blocked  by 
tamoxifen.  Further  changes  were  produced  by  exposure 
and  development  of  resistance  to  fulvestrant  including  a 
loss  of  ERa  Ser118  activation,  which  could  account  for 
its  loss  of  sensitivity  to  estrogen.  These  data  support  the 
view  that  in  the  early  stages  of  resistance,  SERDs  may 


provide  a  useful  therapeutic  option,  but  other 
approaches  will  be  required  when  resistance  has 
developed  to  these  agents. 
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Abstract — Though  supervised  and  unsupervised  analyses  of 
genomic  data  have  been  intensively  studied  in  recent  years,  little 
effort  has  been  made  to  discover  the  structural  information 
contained  in  the  data.  In  this  work,  we  propose  a  stability  analysis 
guided  supervised  clustering  and  visualization  method  aiming  to 
discover  the  hierarchical  structure  in  gene  expression  data,  which 
we  call  the  “tree  of  phenotypes”.  We  applied  the  method  on  two 
multiclass  gene  expression  microarray  data  sets  and  presented 
the  biological  plausibility  of  the  learned  trees.  We  also  tested  the 
multiclass  classifiers  built  on  the  learned  trees  and  demonstrated 
their  good  classification  performance. 

I.  INTRODUCTION 

High  throughput  genomic  data  usually  consists  of  tens  of 
thousands  of  features  and  a  relatively  small  number  of  samples 
of  different  diseases  or  disease  subtypes.  Supervised  analysis 
of  genomic  data  mainly  refers  to  constructing  classification 
schemes  by  learning  from  the  dependencies  between  gene 
expression  values  and  the  given  labels  of  phenotypes  [1]. 
The  common  practice  of  supervised  learning  is  to  design 
and  apply  multiclass  classifiers  with  parallel  structures,  such 
as  multiclass  Support  Vector  Machines  [2],  neural  networks 
[3],  and  Nearest  Shrunken  Centroids  [4].  When  the  number 
of  classes  increases,  the  complexity  of  a  parallel- structured 
classifier  also  increases.  Given  limited  samples  in  a  high 
dimensional  space,  it  is  difficult  to  avoid  overfitting  without 
compromising  the  prediction  performance.  It  is  also  hard  to 
find  a  single  set  of  genes  on  which  the  entire  multiclass 
problem  can  be  reliably  solved. 

Using  tree  classifiers,  we  can  alleviate  these  problems  [5], 
[6],  [7].  In  the  tree  paradigm,  the  original  classification  task 
is  tackled  by  a  series  of  simpler  tasks.  Each  task  might  only 
require  simple  classification  models  to  yield  high  prediction 
accuracy.  The  classification  accuracy  of  tree-based  methods 
is  generally  competitive  with  other  classification  paradigms 
and  structures.  Achieving  good  classification  performance  is 
not  the  only  motivation  for  exploring  the  tree  of  phenotypes. 
The  tree  structure  provides  coarse-to-fine  views  of  the  data 
and  its  class  structure,  which  reveals  the  relationship  between 
classes  (in  our  application  context,  phenotypes);  moreover,  at 
each  node  in  the  tree,  we  can  identify  a  small  set  of  features 
(genes)  that  well  account  for  the  differences  between  the 
phenotypes  associated  with  a  given  node.  The  tree  structure 


can  be  constructed  based  on  a  priori  knowledge,  or  it  can 
be  learned  directly  from  the  data  [5].  For  example,  Shedden 
et  al.  proposed  a  pathological  tree  based  on  tumor  ontogeny 
[6].  Such  a  tree  is  independent  of  any  particular  data  set  and 
therefore  is  not  affected  by  the  sample  size  and  the  data  quality. 
The  main  problem  associated  with  this  approach  is  that  the 
information  expressed  by  the  tree  is  largely  confined  by  the 
prior  knowledge,  which  could  be  incomplete  or  inaccurate. 
For  example  in  Shedden  et  al.’s  tree,  7  tumor  types,  such 
as  lung  cancer  and  prostate  cancer,  are  grouped  together  to 
form  a  single  node  that  represents  the  non-Mullerian  tumors. 
This  does  not  provide  any  information  about  the  relationships 
between  these  seven  types. 

In  this  paper,  we  propose  a  method  that  learns  this  structural 
information  from  gene  expression  profiles  and  their  pheno¬ 
typical  labels.  We  first  introduce  the  stability  analysis  based 
method  for  learning  a  tree  of  phenotypes  from  gene  expression 
data.  We  then  describe  how  to  construct  a  tree  classifier  given 
the  learned  tree  structure.  Lastly,  we  apply  our  method  on  two 
microarray  data  sets,  identify  the  biological  structures  obtained 
by  our  tree  learning  algorithm,  and  demonstrate  the  prediction 
performance  of  the  tree  classifiers. 

II.  LEARNING  THE  TREE  OF  PHENOTYPES 

A  tree  of  phenotypes  is  a  natural  way  to  describe  the 
relationship  between  diseases  or  disease  subtypes.  We  devise 
a  method  that  learns  the  relation  between  phenotypes  with  the 
guidance  of  human  interaction,  namely  Color-Coded  Super¬ 
vised  Mode  Visual  Statistical  Data  Analyzer  (ccsm VISDA), 
an  extension  of  the  original  VISDA  algorithm  [8].  This  method 
gives  the  capability  to  discern  unknown  relationships  between 
phenotypes  that  are  latent  within  the  data.  To  assure  good  gen- 
eralizability  for  small  sample  sizes,  the  tree  learning  procedure 
includes  a  leave-one-out  “stability”  analysis  that  we  propose. 
The  final  predicted  tree  is  the  one  receiving  the  most  votes, 
among  all  the  leave- one- out  trees. 

A.  Learning  Trees  by  ccsmVISDA 

The  ccsmVISDA  algorithm  hierarchically  displays  the 
classes  and  constructs  a  tree.  We  call  a  tree  node  with  two  or 
more  classes  a  composite  node,  and  a  tree  node  with  only  one 
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class  a  terminal  node.  Starting  from  the  root  node,  samples 
are  partitioned  into  clusters  to  grow  the  tree.  A  cluster  is 
considered  as  a  composite  node  if  it  contains  more  than  one 
class;  otherwise,  it  is  a  terminal  node.  At  each  composite  node, 
the  local  data  is  first  projected  onto  a  visualization  subspace 
that  allows  the  user  to  interactively  initialize  the  clustering. 
The  cluster  partition  is  iteratively  updated  until  a  stable  state 
is  reached.  During  the  updating,  samples  from  the  same  class 
are  forced  to  be  assigned  to  the  same  cluster,  i.e.  clusters  learn 
to  fully  “own”  either  one  or  multiple  classes. 

Before  constructing  the  tree,  we  first  filter  the  genes  by  their 
signal-to-noise  ratios  (SNR).  The  purpose  is  to  remove  those 
non-discriminatory  genes  and  ease  the  computational  demand. 
Suppose  the  data  set  consists  of  K  classes  with  p  genes;  each 
class  has  samples,  k  =  1, . . .  ,  FT.  Denote  the  mean  and 
standard  deviation  of  data  from  class  k  and  gene  i  by  /i^  and 
<Jik ,  where  k  =  l,...,K  and  %  =  1, . .  *  ,p.  We  define  the  SNR 
via: 


SNR(i)  = 


Em=1  '^2y=u+l7ru'^v{^iu  k'iv)2 


EK  2 

k=l  'KkO’ik 

nk  =  nk/Ylf=tnji  k  =  l,...,K. 


(1) 


The  top  m  genes,  with  highest  SNRs,  are  used  to  represent  the 
data.  Here  m  is  proportional  to  K  and  determined  empirically. 
We  apply  ccsmVISDA  on  the  filtered  data. 

Suppose  at  a  composite  node  there  are  no  samples  with  m 
genes  coming  from  K0  classes.  Denote  the  mean  vector  and 
the  covariance  matrix  of  class  c  by  fic  and  Cc.  All  the  samples 
are  first  projected  onto  a  two-dimensional  space  selected  by 
multiclass  Fisher’s  discriminant  analysis  [9],  which  utilizes  the 
class  information  to  find  the  most  discriminatory  subspace  for 
the  Kq  classes.  The  projection  space  is  spanned  by  the  two 
vectors  that  maximize  Fisher’s  criterion  [9],  i.e. 


K0)  and  the  initial  location  p^k  for  the  center  of  each  cluster 
k  in  the  projection  plot.  Each  cluster  is  modeled  by  a  single 
Gaussian  distribution.  Denote  the  probability  density  function 
of  a  Gaussian  distribution  by 

p(x|/x,  C),  (7) 


pi  the  mean  vector  and  C  the  covariance  matrix.  In  order  to 
get  a  more  robust  partition  of  the  samples,  the  user  is  further 
required  to  select  two  more  partition  schemes  that  have  M  —  1 
and  M  +  1  clusters.  All  three  partition  schemes  will  undergo 
the  same  clustering  process.  For  the  partition  scheme  with  M0 
clusters  (M0  G  {M  —  1,  M,  M+l}),  after  the  user  has  selected 
the  centers,  each  sample  x*  is  assigned  to  a  cluster  gi  such 
that 

9i  =*  argmin  {||x*  -  Mxfc||}  •  (8) 


This  initial  partition  of  the  data  into  clusters  is  iteratively 
updated  by  an  EM-like,  two  step  procedure.  Denote  the 
partition  at  the  n’th  iteration  by 

S(n>  =  {5|n),5<n),...,5^o)},  (9) 

where  is  the  index  set  of  the  samples  that  are  assigned 
to  cluster  k. 

In  the  first  step  of  each  iteration,  each  sample  is  assigned 
to  a  cluster.  Define 


z 


(n) 

ik 


5{gi  -  k), 


r(n) 


in)  n(n)s 
FxA  ?  ^xfc  ) 


|M: 


(n)  p(n)\ 


n  =  0 
n  >  1 


where  S(-)  is  the  Kronecker  delta  function,  i.e. 


(10) 


<%)  = 


1,  0  =  0 

0,  0/0. 


Wo  =  argmax  {trace(WTSw1SbW)}  ,  (2) 

w 

where  Wo  is  a  m  by  2  matrix.  Here,  the  within  class  scatter 
matrix  Sw  is  defined  as  [9] 

K0 

S\V  =  ^7TcCc  (3) 

c=  1 

and  the  between  class  scatter  matrix  is  defined  as  [9] 

K0-l  K0 

Sb  =  ^  ^  ^  ^  KiTTj  (y  f-li  /i-j)(/i^  Mj)  •  (4) 

7=1  j=i+ 1 

Here  nc  is  the  mixing  proportion  of  class  c,  i.e. 

7vc  =  \Ic\/no  (5) 

with  Ic  the  index  set  of  the  samples  from  class  c,  and  \IC\  the 
size  of  set  Ic.  Each  sample  t  in  the  m-dimensional  space  is 
projected  into  the  2-D  space 

x  =  Wjt.  (6) 

Given  the  2-D  visualization  of  the  samples,  the  user  is 
required  to  determine  both  the  number  of  clusters  (M,  M  < 


(n) 

When  n  >  1,  zik  is  the  a  posteriori  probability  that  the 
sample  x*  belongs  to  cluster  k.  Each  class  l  is  assigned  to  a 
cluster  ki  such  that 


h 


arg  max 


(11) 


Note  that,  through  this  operation,  multiple  classes  may  be 
assigned  to  the  same  cluster.  The  partition  S ^  is  updated 
accordingly: 

5f)=  U  7i’  k  =  l,...,M0.  (12) 


In  the  second  step,  the  mean  vector  and  covariance  matrix 
of  each  cluster  are  updated  by 


k::i]  = 


^ i£S 


(n)  Xi 


?(n) 


(n+1) 


M0  Q(n) 


Ei=°i  \s: 


(13) 


(14) 
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EieS(n)(Xi  -  -  Mifc+1))T 

- Fn - 


(15) 


of  the  output  tree  structures  associated  with  the  tree  learning 
algorithm.  Let 

T  :  Z+,  (18) 


Eqs.  10-15  describe  the  operations  for  updating  the  partition. 
They  are  repeated  until  the  partition  (Eq.9)  no  longer  changes. 

The  best  partition  is  then  determined  by  the  user,  with  the 
guidance  of  the  Minimum  Description  Length  (MDL)  score 
for  each  of  the  three  partitions  [8].  The  MDL  score  for  the 
partition  with  M0  clusters  is 

no  j\/r  -| 

MDL(M0)  =  -  ^  logp(xj)  H - | - logn0,  (16) 

1=1 

M0  e  {M  —  1,M,  M  +  1},  where  p(x)  is  the  probability 
density  function  of  the  mixture  of  Gaussians 

M0 

p(x)  =^7Tk  p(x  |  /Xxfc,  Cxjfc).  (17) 

k=  1 

After  the  optimal  partition  is  determined,  each  cluster  be¬ 
comes  a  child  node  of  the  current  composite  node,  and  the 
partitioning  is  repeated  until  each  node  contains  all  the  samples 
from  one  class. 

We  pursue  such  a  human-interactive  visualization  approach 
within  the  mixture  modeling  for  two  complementary  reasons. 
First,  users  of  ccsm-YISDA  can  incorporate  their  knowledge 
about  the  relations  between  phenotypes  during  construction  of 
the  trees;  second,  such  an  approach  can  help  users  to  discover 
new  relations  between  phenotypes.  The  merits  of  this  approach 
have  been  empirically  evidenced  by  biological  studies  [10], 
[11]. 

B.  The  Stable  Solution  of  the  Tree 

Learning  the  tree  from  all  the  samples  may  cause  overfitting 
due  to  small  sample  size  and  potentially  poor  data  quality. 
Thus,  we  embed  the  ccsmVISDA  learning  procedure  within 
a  leave-one-out  stability  analysis  to  generate  “leave-one-out 
trees”,  and  then  take,  as  the  final  solution,  the  tree  in  this  set 
whose  hierarchical  class  structure  has  the  highest  frequency 
of  occurrence.  This  winning  tree  reflects  the  underlying  stable 
structural  information  in  the  data  because  it  is  learned  from 
the  data  set  and  amongst  all  learned  trees  best  survives  small 
disturbances  of  the  data  [12].  The  stability  analysis  based 
ccsmVISDA  algorithm  (SA-ccsmVISDA)  is  more  robust  than 
other  solutions  in  the  sense  that,  given  different  realizations 
of  the  data  distribution,  the  algorithm  will  output  similar 
solutions.  The  robustness  of  the  tree  learning  algorithm  is 
critical  to  making  scientific  discoveries,  since  the  learned  tree 
must  be  reproducible  with  high  probability  in  the  face  of  small 
data  variations,  in  order  to  be  used  as  a  hypothesis  for  the 
underlying  relations  between  phenotypes. 

Let  the  set  of  all  possible  hierarchical  class  structures  with 
K  terminal  nodes  be  the  sample  space  0#,  which  contains  all 
distinguishable  outputs  of  a  tree  learning  algorithm  on  a  given 
set  of  samples  consisting  of  K  classes.  There  is  a  distribution 


where  T  is  a  random  variable  whose  values  are  the  indices  of 
tree  structures.  We  can  use  entropy  to  measure  the  stability  of 
the  tree  learning  algorithm,  given  by: 

H(L)  =  -  D  p(t)\ogp{t),  (19) 

tez+ 

where  p(t)  is  the  probability  mass  function  of  the  derived 
distribution  for  random  variable  T;  L  is  the  tree  learning  al¬ 
gorithm.  We  can  estimate  the  distribution  p(t)  by  the  empirical 
distribution  of  the  leave-one-out  tree  structures  generated  by 
SA-ccsmVISDA.  Without  loss  of  generality,  we  can  always 
define  T  such  that  it  maps  the  distinctive  outputs  of  SA- 
ccsmVISDA  to  a  set  of  consecutive  positive  integers  starting 
from  1.  We  simply  set  p(t)  equal  to  0  for  those  structures  that 
do  not  appear  in  the  leave-one-out  trees.  In  the  experiments 
we  plot  such  empirical  distributions  and  calculate  the  stability 
measure  using  Eq.19. 

In  the  experiments,  we  also  show  the  distance  measures 
along  the  winning  tree  structure.  The  distance  is  defined  as  the 
average  of  the  Fisher’s  criteria  between  all  pairs  of  classes  at  a 
composite  node  in  the  2-D  projection  space,  and  is  calculated 
using  all  the  samples.  The  distance  is  represented  by  the 
common  length  of  the  links  between  a  composite  node  and 
all  its  child  nodes. 

The  stability  analysis  based  ccsm-VISDA  approach  for  tree 
learning  based  on  leave-one-out  is  practically  feasible  given 
the  small  sample  sizes  in  most  of  the  existing  microarray 
data  sets.  The  amount  of  human  interaction  required  for  each 
of  the  leave-one-out  trials  is  not  always  the  same.  In  the 
experiments  we  have  observed  that  leaving  one  sample  out 
each  time  usually  will  not  change  the  structure  of  trees  at 
the  top  levels.  The  projections  will  mostly  change  at  the 
deeper  levels,  in  the  composite  nodes  that  include  the  class 
whose  samples  are  left  out.  For  such  instances,  we  only  need 
to  apply  ccsm-VISDA  once  on  those  branches  that  are  not 
subject  to  change  in  a  specific  set  of  leave-one-out  trials.  Thus, 
in  practice,  the  human  interaction  in  stability  analysis  based 
ccsm-VISDA  is  less  intensive  than  that  required  for  repetition 
of  the  ccsm-VISDA  procedure  over  all  leave-one-out  training 
sets.  If  the  sample  sizes  are  increased  in  future  microarray 
data  sets,  we  can  modify  our  approach  to  be  semi  or  fully 
automated  by  exhaustively  searching  the  optimal  number  of 
sub-clusters  at  each  node,  but  in  this  case  human  users  will 
have  less  or  no  chance  to  incorporate  their  knowledge  into  the 
structure  of  the  trees.  For  existing  microarray  data  with  small 
sample  sizes,  our  stability  analysis  based  ccsm-VISDA  method 
gives  a  reasonable  balance  between  robustness  and  practical 
feasibility. 

III.  TREE  CLASSIFIER 

The  tree  structure  learned  by  SA-ccsmVISDA  can  be  used 
to  build  tree-based  classifiers.  Classification  trees  are  a  general 
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framework  for  solving  multiclass  classification  problems  [5]. 
We  can  use  potentially  any  classifier  as  the  node  classifier,  and 
use  a  different  subspace  for  the  particular  classification  task 
on  each  node.  Subspace  feature  selection  will  help  not  only 
improving  classification  performance,  but  also  in  finding  the 
genes  that  primarily  account  for  the  similarities  or  differences 
between  subsets  of  phenotypes  [6]. 

In  the  experiments,  we  used  the  feature  filtering  and  selec¬ 
tion  method  proposed  by  Shedden  et  al.  [6].  First  the  control 
genes  are  removed  from  the  data.  Then  the  expression  values 
are  transformed  by  log[max(x,  0)  +  50].  All  those  genes  for 
which  the  standard  deviation  across  all  samples  is  less  than  0.7 
are  removed.  For  each  sub-cluster  on  a  composite  node,  the 
set  of  a  genes  that  have  the  highest  mean  expression  values 
(and  greater  than  the  mean  expression  values  of  all  the  other 
sub-clusters  combined)  are  selected.  All  these  sets  of  genes 
are  combined  to  form  the  subspace.  Here  the  range  of  a  is 
determined  by  the  user  and  all  the  nodes  use  the  same  value 
of  a.  The  optimal  value  of  a  is  selected  to  achieve  the  lowest 
leave-one-out  cross  validation  error  rate  of  the  tree  classifier. 

We  use  one-versus-rest  multiclass  Support  Vector  Machines 
(OVR-MSVM)  [2]  as  the  node  classifiers.  For  a  given  data 
set,  we  evaluate  both  the  hard  classifier  whose  outputs  are 
class  labels  and  the  soft  classifier  whose  outputs  are  confidence 
values  [13]. 

In  the  hard  classification  scheme,  for  each  sub-cluster  on 
a  composite  node  a  binary  SVM  is  trained  in  the  selected 
subspace  to  separate  the  sub-cluster  from  all  others.  A  testing 
sample  will  be  assigned  to  the  sub-cluster  whose  associated 
binary  SVM  has  the  largest  real- valued  output,  and  will  be 
passed  to  the  child  node  corresponding  to  the  sub-cluster. 
When  the  testing  sample  reaches  a  terminal  node,  it  will  be 
classified  to  the  phenotype  associated  with  the  terminal  node. 

In  the  soft  classification  scheme,  an  OVR-MSVM  is  trained 
in  the  same  way  as  in  the  hard  classification  scheme  except 
that  for  each  binary  SVM  the  real-valued  output  is  transformed 
to  produce  a  posteriori  probabilities.  We  applied  the  method 
proposed  by  Platt  [13]  to  derive  this  transformation.  The  a 
posteriori  probability  output  of  an  OVR-MSVM  is  given  in 
the  form 

=  Mo(X)  *  =  1, .  •  • ,  Mq,  (20) 
Efc=i<?fc(x) 

where  <^(x)  is  the  transformed  SVM  output  specified  in  [13]. 
When  a  testing  sample  is  tested  from  the  root  node  to  each 
terminal  node,  the  simulated  a  posteriori  probabilities  are 
multiplied  together.  The  output  at  each  terminal  node  is  taken 
as  the  a  posteriori  probability  of  the  sample  belonging  to  the 
phenotype  associated  with  the  node.  The  sample  is  assigned 
to  the  phenotype  with  the  highest  simulated  a  posteriori 
probability.  We  will  see  in  the  experiments  that  the  soft 
classification  scheme  improves  the  performance  of  the  tree 
classifier. 

IV.  EXPERIMENTS 

In  the  experiments,  we  applied  the  stability  analysis  based 
ccsmVISDA  on  the  Muscular  Dystrophy  data  set  to  generate 


the  tree  of  9  subtypes  of  muscular  dystrophies  and  on  the  MIT 
cancer  data  to  generate  the  tree  of  14  cancers.  For  each  data 
set,  we  also  evaluated  the  performance  of  the  tree  classifiers 
built  on  the  tree  of  phenotypes. 

A.  Muscular  Dystrophy  Data  Set 

The  muscular  dystrophy  data  set  (provided  by  Children  Na¬ 
tional  Medical  Center  (CNMC),  Center  for  Genetic  Medicine) 
consists  of  108  samples  with  11252  genes  from  9  diagnostic 
groups  of  muscular  dystrophies.  The  name  and  the  number  of 
samples  of  each  group  are:  amyotrophic  lateral  sclerosis  (ALS, 
n=9);  acute  quadriplegic  myopathy  (AQM,  n=5);  calpain  III 
deficiency  (Calpain3,  n=10);  Duchenne  muscular  dystrophy 
(DMD,  n=10);  dysferlin  deficiency  (Dysferlin,  n=10);  fukutin 
related  protein  deficiency  (FKRP,  n=7);  fascio scapulohumeral 
dystrophy  (FSHD,  n=14);  normal  human  muscle  (NHM, 
n=18);  and  juvenile  dermatomyositis  (JDM,  n=25). 

We  applied  SA-ccsmVISDA  on  the  data  and  derived  108 
leave-one-out  trees  showing  12  different  structures.  The  em¬ 
pirical  distribution  of  the  trees  is  shown  in  Fig.l.  The  entropy 
calculated  using  Eq.19  is  about  1.4208.  The  frequency  of  the 
winning  tree  is  67/108  ~  0.62. 


Tree  ID 

Fig.  1.  The  empirical  distribution  of  the  tree  structures  for  the  Muscular 
Dystrophy  data  set.  The  abscissa  is  the  index  of  the  structures,  which  is  in 
descending  order  of  frequencies.  The  ordinate  is  the  frequency.  The  number 
on  the  top  of  each  bar  is  the  number  of  occurrences  of  the  structure.  The 
entropy  H(L)  of  the  set  of  tree  structures  is  about  1.4208. 

The  structure  of  the  winning  tree  is  shown  in  Fig. 2  with 
subtype  names  on  the  terminal  nodes.  The  winning  tree  is 
supported  by  many  known  clinical,  genetic  and  histological 
features  of  these  disorders.  ALS  is  the  only  denervating 
disorder,  due  to  die-back  of  motor  neurons.  A  number  of 
the  muscular  dystrophies  are  caused  by  abnormalities  in 
the  plasma  membrane  of  the  muscle  fiber:  Calpain3,  DMD, 
Dysferlin  and  FKRP  are  all  such  membrane  dystrophies,  and 
all  group  together.  FSHD  is  a  unique  disorder  due  to  a 
heterozygous  deletion  in  chromosome  4q,  and  by  our  SA- 
ccsmVISDA  approach  this  maps  distinctly,  as  does  normal 
human  muscle  (NHM)  and  an  autoimmune  disease,  JDM. 

Based  on  the  winning  tree  structure,  we  built  hard  and 
soft  tree  classifiers  using  the  subspace  selection  method  [6] 
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Fig.  2.  The  structure  of  the  winning  tree  for  the  Muscular  Dystrophy  data 
set.  The  subtypes  of  muscular  dystrophies  are  shown  on  the  terminal  nodes. 
The  distance  measures  are  shown  along  the  links. 

TABLE  I 

Leave-one-out  error  rates  (in  percentage)  for  hard  and  soft 

CLASSIFICATION  SCHEMES  ON  THE  MUSCULAR  DYSTROPHY  DATA.  IN 
PARENTHESES  ARE  THE  VALUES  OF  a.  THE  LOWEST  ERROR  IS  IN  BOLD 
FACE. 


Fig.  3.  Leave-one-out  error  rates  of  hard  and  soft  classifiers  on  the  Muscular 
Dystrophy  data  set  for  C  =  1.0. 


C 

Hard 

Soft 


i  o.ooi  o.oi _ or _ l.o  lo.o 

23.15  (800)  12.96  (250)  12.04  (50)  12.04  (50)  12.04  (50) 
12.04  (300)  11.11  (400)  11.11  (300)  10.19  (400)  11.11  (50) 


described  in  section  3  and  OVR-MSYM  as  the  node  classifier. 
The  gene  expression  values  were  transformed  by  log(x)  before 
evaluating  the  classifiers.  The  OVR-MSVM  uses  linear  SVMs 
as  the  binary  components.  The  complexity  of  a  linear  SVM 
can  be  controlled  by  a  penalty  value,  C.  We  use  the  same 
C  value  for  each  linear  SVM  in  each  of  the  OVR-MSVMs. 
We  tested  5  values  for  C\  0.001,  0.01,  0.1,  1.0  and  10.0. 
In  order  to  determine  an  optimal  subspace  size,  for  each  C 
value  we  tested  25  different  values  for  a:  1,2,  5,  10,  20  and 
the  sequence  {50 k,k  =  1, . . . ,  20}.  All  these  tree  classifiers 
are  evaluated  by  leave-one-out  cross  validation.  In  Tabled,  we 
list  the  lowest  error  rates  (in  percentage)  for  both  soft  and 
hard  classification  schemes  for  each  C  value.  The  values  in 
parentheses  are  the  values  of  a  at  which  the  performances  are 
achieved.  In  Fig. 3,  we  plot  the  leave-one-out  error  rates  of  soft 
and  hard  classification  as  functions  of  a  for  C  =  1.0.  It  can 
be  seen  that  soft  classification  improves  the  performance  for 
most  cases.  The  lowest  error  rate  is  10.19%  when  C  =  1.0 
and  a  =  400. 

B.  MIT  Cancer  Data  Set 

The  MIT  cancer  data  was  originally  proposed  in  Ra- 
maswamy  et  al.  [2].  The  data  set  has  16063  genes.  It  contains 
a  training  set  with  144  samples  and  an  independent  testing  set 
with  54  samples  of  14  cancer  types.  The  control  genes  are 
removed  from  the  data  before  all  the  following  experimental 
steps. 

We  applied  SA-ccsmVISDA  on  the  training  set  and  gener¬ 
ated  144  trees  in  the  leave-one-out  loop.  The  144  trees  demon¬ 
strate  20  different  structures,  whose  empirical  distribution  is 
shown  in  Fig.  4.  The  entropy  of  the  empirical  distribution 
calculated  using  Eq.19  is  about  1.3344.  The  frequency  of  the 
winning  tree  is  102/144  ~  0.71.  The  structure  of  the  winning 
tree  with  the  cancer  types  shown  on  the  terminal  nodes  is 


illustrated  in  Fig.  5. 


1  2  3  4  5  6  7  8  9  10  1 1  12  13  14  15  16  17  18  19  20 
Tree  ID 


Fig.  4.  The  empirical  distribution  of  the  tree  structures  for  the  MIT  cancer 
data.  The  entropy  H(L)  of  the  set  of  tree  structures  is  about  1.3344. 

Using  the  same  scheme  as  for  the  muscular  dystrophy  data 
set,  we  evaluated  the  hard  and  soft  tree  classifiers  for  the  MIT 
cancer  data.  Before  the  evaluation,  the  gene  expression  values 
were  transformed  and  the  genes  were  filtered  just  as  Shedden 
et  al.  did  in  [6].  The  lowest  leave-one-out  error  rates  for 
different  values  of  C  are  listed  in  Table.II.  The  curves  of  leave- 
one-out  errors  of  soft  and  hard  classifications  as  functions  of 
a  are  illustrated  in  Fig.  6.  The  best  classification  accuracy  of 
our  tree  classifiers  is  about  87.5%  when  C  =  1  and  a  =  200. 
Our  results  compare  favorably  with  those  of  Ramaswamy  et 
al.  [2],  who  used  a  parallel  OVR-MSVM  to  classify  the  14 
cancers  and  achieved  a  78%  accuracy  on  the  training  set  of 
144  samples  using  leave-one-out  cross  validation. 

As  a  means  to  explore  the  biological  implications  of  our 
solution,  we  compared  our  tree  to  Shedden  et  al.’s  tree  based 
on  pathologic  and  ontologic  knowledge  [6].  Our  solution 
has  notable  similarities  to  this  tree,  consistently  classifying 
lymphoma,  leukemia,  CNS  and  epithelial  cancers  into  groups 
in  which  lymphoma  and  leukemia  are  closely  related  and 
CNS  and  epithelial  cancers  are  closely  positioned.  Cancers  of 
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Fig.  5.  The  structure  of  the  winning  tree  for  the  MIT  cancer  data.  The 
cancer  types  are  shown  on  the  terminal  nodes.  The  distance  measures  are 
shown  along  the  links. 
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Fig.  6.  Leave-one-out  error  rates  of  the  hard  and  soft  classifiers  on  MIT 
cancer  data  for  C  =  1.0. 


Leave-one-out  error  rates  (in  percentage)  for  soft  and  hard 

CLASSIFICATION  SCHEMES  ON  THE  MIT  CANCER  DATA.  IN  THE 
PARENTHESES  ARE  THE  VALUES  OF  a.  THE  LOWEST  ERROR  IS  IN  BOLD 
FACE. 


within  the  stability  analysis  framework,  we  can  generate  more 
stable  and  more  generalizable  tree  solutions  given  a  limited 
number  of  samples. 


c 

0.001  0.01  0.1  1.0  10.0 

Hard 

Soft 

21.53  (450)  16.67  (200)  15.97  (200)  15.97  (200)  15.97  (200) 
16.67  (250)  13.19  (250)  13.19  (200)12.50  (200)  13.19  (200) 

the  uterus,  breast,  lung,  colon,  bladder,  kidney  and  pancreas 
are  consistently  and  appropriately  classified  into  the  same 
group.  A  major  difference  between  the  two  trees  is  the 
location  of  the  mesotheliomas  (a  rare  cancer)  and  melanomas. 
Applying  an  independent  data-driven  approach  to  the  same 
data  set,  Tibshirani  and  Hastie  [7]  generated  a  tree  broadly 
similar  to  ours,  in  which  the  melanomas  and  mesotheliomas 
also  clustered  differently  than  predicted  by  the  Shedden  et 
al.  construct.  Further  analysis  of  the  similarities  among  the 
melanomas,  mesotheliomas  and  their  most  closely  related 
cancers  may  generate  new  insights  into  common  molecular 
functions  among  these  cancers. 

Y.  CONCLUSION 

In  this  paper,  we  proposed  a  stability  analysis  based  data 
visualization  algorithm  that  learns  the  tree  of  phenotypes  from 
genomic  data  with  human  guidance.  We  applied  the  algorithm 
on  two  gene  expression  microarray  data  sets.  The  efficacy  of 
the  algorithm  is  demonstrated  by  the  biological  information 
discerned  in  the  derived  trees.  We  also  demonstrated  the 
prediction  performance  of  the  multiclass  classifiers  built  on 
the  derived  trees. 

This  hierarchical  representation  of  phenotypes  has  the  power 
to  reveal  both  global  and  local  structures  that  are  important 
for  understanding  the  relationships  between  phenotypes.  By 
selecting  a  different  subspace  on  each  composite  node,  we  can 
find  the  genes  that  are  important  for  explaining  either  the  sim¬ 
ilarity  or  the  difference  between  phenotypes  in  a  given  group. 
The  human  guided  visualization  approach  is  more  robust  to 
corrupted  data  and  the  small  sample  size  problem  than  purely 
automated  methods.  By  embedding  the  tree  learning  procedure 
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Abstract  -  In  this  paper,  we  report  a  new  gene  clustering 
approach  -  non-negative  independent  component  analysis 
(nICA)  -  for  microarray  data  analysis.  Due  to  positive  nature  of 
molecular  expressions,  nICA  fits  better  to  the  reality  of 
corresponding  putative  biological  processes.  In  conjunction  with 
nICA  model,  Visual  Statistical  Data  Analyzer  (VISDA)  is 
applied  to  group  genes  into  modules  in  the  latent  variable  space. 
The  experimental  results  show  that  significant  enrichment  of 
gene  annotations  within  clusters  can  be  obtained. 

Keywords  -  non-negative  ICA,  latent  variable  model,  gene 
clustering,  module  discovery,  microarray  data  analysis 

L  Introduction 

Microarray  technologies  provide  powerful  tools  for 
genome- wide  measurement  of  gene  expressions.  To 
discover  functional  modules  involved  in  pathway  signaling  or 
gene  regulation,  new  computational  methods  are  required  for 
modeling  and  analysis  of  microarray  data  of  interest  [1]. 

Gene  clustering  is  widely  used  in  the  analysis  of  gene 
expression  data  by  partitioning  genes  into  clusters  sharing 
similar  expression  patterns.  The  underlying  assumption  is  that 
genes  with  similar  patterns  are  more  likely  associated  with 
common  functions.  Hierarchical  clustering  and  self¬ 
organizing  maps  [2],  have  been  applied  to  group  the  genes 
into  functional  modules.  Recently,  Independent  Component 
Analysis  (ICA)  has  been  proposed  for  modeling  gene  clusters 
[3].  In  contrast  to  traditional  clustering  methods,  ICA-based 
clustering  relies  on  a  linear  combination  of  latent  biological 
processes  and  has  revealed  the  gene  clusters  with  significant 
enrichment  of  gene  annotations  or  functional  categories  [3]. 
In  contrast  to  PCA,  ICA  decomposes  input  data  into 
components  as  independent  as  possible,  showing  some 
advantages  over  PCA  for  gene  module  decomposition  [5]. 

In  this  paper,  we  report  the  application  of  non-negative 
ICA  (nICA)  for  gene  clustering,  exploiting  the  non-negative 
nature  of  molecular  expressions.  In  principle,  nICA  can  be 
thought  as  a  projection  method  where  the  expression  levels 
are  projected  onto  some  new  non-negative  bases  (i.e. 
components)  with  minimum  statistical  dependence.  The  nICA 
representation  shall  better  reflect  the  biological  reality.  We 
then  use  Visual  Statistical  Data  Analyzer  (VISDA)  [6]  to 
generate  gene  modules  in  the  latent  variable  space.  VISDA 
uses  hierarchical  Standard  Finite  Normal  Mixtures  (SFNM) 
to  model  clustered  data  where  each  gene  belongs  to  each 
cluster  with  a  posterior  probability.  The  clustering  procedure 
follows  a  hierarchy  fasion.  At  each  level  of  the  hierarchy, 
each  cluster  is  considered  for  further  split,  until  no  cluster  is 
decomposable  according  to  the  Minimum  Description  Length 
(MDL)  criterion  or  human  justification. 


This  paper  is  organized  as  follows.  In  section  II,  we 
introduce  the  principle  of  nICA  for  finding  gene  module 
composites  and  a  gradient  descent  algorithm  of  nICA.  A  brief 
description  of  the  VISDA  algorithm  is  also  given  in  Section 
II  to  cluster  independent  components  as  a  post-processing  of 
nICA  modeling.  The  application  of  nICA  and  VISDA  to 
yeast  data  will  be  reported  in  Section  III.  Discussions  and 
conclusions  are  given  in  Section  IV. 

II.  Methodology 

The  problem  of  basic  ICA  is  given  according  to  the 
following  linear  relation: 

x  =  As  (1) 

where  s  =  (si,  s2...  sN)  is  a  vector  of  real  independent  sources 
and  x  =  (x1?  x2...  xN)  is  an  observation  vector.  The 
assumptions  of  ICA  are  that  sources  are  mutually 
independent  and  of  non-Gaussian  distribution  except  for  at 
most  one  source.  When  we  apply  ICA  to  real-world  problems 
like  gene  expression  analysis,  the  situation  is  different  from 
the  above  ideal  case  because  of  the  ambiguities  of  ICA:  the 
sign  and  permutation  of  sources.  To  resolve  these 
ambiguities,  many  researchers  make  further  assumptions  to 
constrain  the  ICA  model.  For  example,  non-negativity  is  a 
natural  constraint  for  many  real-world  applications,  such  as 
blind  separation  of  natural  scenes  [7].  Since  we  assume  that 
the  underlying  biological  processes  are  independent  and  their 
expression  levels  should  be  non-negative,  nICA  is  believed  to 
be  a  more  proper  model  to  represent  a  linear  influence  of 
hidden  cellular  variables  than  ICA  is.  By  projecting  the  data 
to  the  latent  space  spanned  by  these  non-negative 
independent  processes,  fine  structure  of  co-regulation  of 
genes  is  maintained  and  made  prominent.  VISDA  clustering 
is  then  applied  to  catch  the  characteristics  of  those  subtle 
differences,  which  may  lead  to  identify  more  coherent  gene 
groups  (Fig.  1). 

A.  nICA-based  decomposition  and  the  algorithm  of  nICA 

As  it  has  been  known,  clustering  by  expression  pattern  or 
“co-expressed”  genes  under  limited  experimental  conditions 
does  not  provide  the  best  possible  grouping  of  genes  by 
biological  processes  [8].  ICA-based  gene  clustering  approach, 
on  the  other  hand,  is  built  upon  a  latent  variable  model  of 
gene  module  composite.  The  attraction  of  ICA  clustering  lies 
that  it  can  account  for  independent  hidden  effects  that 
influence  gene  expression.  When  we  introduce  the  non¬ 
negativity  into  the  ICA  algorithm,  the  resulting  nICA 
approach  can  incorporate  prior  knowledge  for  better 
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modeling  hidden  sources  while  keeping  all  the  advantages  of 
ICA  approach. 

In  our  nICA  model,  gene  express  is  a  linear  combination 
of  biological  mechanisms  including  pathways  of  signaling 
substances,  transcription  factors  and  their  binding  sites  in  the 
promoter  regions  of  genes,  as  well  as  other  different  kinds  of 
regulation  [3].  We  use  nICA  to  project  expression  data  X  to 
the  independent  mode  in  order  to  highlight  these  factors.  We 
assume  that  x(i,  j)  which  is  the  expression  level  of  gene  i 
under  phenotype  j  is  expressed  by  the  sum  of  non-negative 
independent  putative  biological  processes  Sk(i),  k  =  1,  2,  ..., 
N,  weighted  by  the  involvement  strength  ak(j),  k  =  1,  2,  . . .,  N. 


^  Pre-screening 


Fig.  1  Framework  of  nICA  and  VISDA  for  the  composite  module  discovery 

In  [7],  the  author  suggested  a  mathematical  approach  to 
impose  non-negative  constraint  on  sources.  And  if  we  define 
non-negative  well-grounded  sources  as: 

p(sk  <  5)  >  0  for  \/  8  >  0 

p(sk  <  0)  =  0  k  =  1,2,...,  A  (2) 

then  it  has  been  proven  [7]  that  we  can  find  y  =  Us,  where  U 
is  a  square  orthonormal  rotation  and  permutation  matrix,  i.e. 
the  elements  of  y-x  of  y  are  a  permutation  of  sources  if  and 
only  if  all  yis  are  all  non-negative.  We  note  that  y  =  Us  can  be 
re-written  as  y  =  Wz  =  WVx  =  WVAs  with  V  a  whitening 
matrix,  z  the  pre-whitened  observation  vector  and  W  an 
unknown  orthonomal  (rotation)  matrix.  Therefore  we  can 
consider  nICA  as  a  procedure  with  the  following  two  steps: 

1)  remove  the  second  order  statistics  by  whitening;  2)  search 
for  a  rotation  matrix  where  all  the  data  fit  into  the  positive 
quadrant. 

As  described  in  [7],  we  can  use  the  cost  function  J  defined 
in  the  following  to  find  the  global  minimum: 

J(W)  =  £{|z-WV||2}  (3) 

y  =  Wz 

y*  =  max(0,y.) 

y+=(yfyfL  y+N) 

Based  on  the  gradient  descent  rule,  a  learning  algorithm  to 
find  the  de-mixing  matrix  W  is  defined  as  follows  [7] : 

1)  Pre-whitening  the  observed  data  x: 

z  =  Vx  (4) 

2)  Using  gradient  descent  algorithm  to  minimize  the  cost 
function  (3): 

w  =  W  (5) 


3)  Projecting  the  unconstraint  gradient  descent  set  onto  a  set 
of  orthonormal  vectors. 

B.  Pre-screening  for  the  clustering 

After  nICA,  we  obtain  some  independent  components 
describing  some  distinct  biological  processes.  In  these 
putative  biological  processes,  some  genes  showing  relatively 
high  or  low  expression  levels  are  most  interesting.  We  will 
use  a  pre-screening  procedure  to  single  out  these  genes. 

Specifically,  we  can  select  a  subsets  of  genes  within  one 
of  the  components,  which  includes  over  expressed  genes 
(which  are  activated)  and  down  expressed  genes  (which  are 
repressed)  according  to  the  value  of  each  gene  in  the 
component  [3]: 

a  subset  of  genes  = 

{genes  \fgenes  e  C%  of  largest  values  ofyj 

U {genes  \fgenes  e  C%  of  smallest  values  of y.} 

By  this  pre-screening  step,  we  actually  remove  some 
invariant  genes  in  each  component.  By  taking  the  union  of  the 
selected  genes,  we  provide  a  pool  of  more  meaningful  and 
relevant  genes  to  biological  processes  for  the  next  step¬ 
clustering  -  to  identify  genes  that  belong  to  co-expressed 
modules  in  each  component. 

C.  VISDA  clustering 

In  this  step,  we  will  cluster  genes  into  modules  associated 
with  their  values  in  the  independent  components.  VISDA 
employs  the  hierarchical  SFNM  model  for  hierarchical 
clustering.  The  hierarchical  SFNM  model  uses  the  following 
probability  density  function  to  describe  the  relationship 
between  successive  levels  in  the  hierarchy, 

f(rM)=IXI>#g(rle#) 

4>  4  v  ’ 

IX =i  and  I>#=1 

k= 1  7=1 

where  r  denotes  the  genes  to  be  grouped,  the  upper  level  has 
K0  clusters,  the  Ml  cluster  in  the  upper  level  has  LK  sub¬ 
clusters  in  the  lower  level,  7ik  is  the  mixing  proportion  of  the 
Ml  up  level  cluster,  n )\k  is  the  mixing  proportion  of  the  yth 
sub-cluster  in  the  Ml  upper  level  cluster,  g(«)  is  Gaussian 
distribution  function,  0^  are  the  parameters  associated  with 
the  sub-cluster.  The  fitting  process  of  this  model  is  executed 
by  the  Expectation  Maximization  (EM)  algorithm  [6],  which 
achieves  a  local  maximum  of  the  likelihood  function. 

For  each  cluster  at  a  level  of  the  hierarchy,  VISDA  uses 
two  different  projection  methods,  Principle  Component 
Analysis  (PCA)  and  Principle  Component  Analysis  - 
Projection  Pursuit  (PCA  -  PPM)  [6],  to  visualize  the  sub¬ 
clusters  within  the  clusters.  The  user  chooses  one  of  the 
projections  that  he/she  thinks  better  revealing  the  data 
structure.  On  the  chosen  projection,  user  initializes  models 
with  different  number  of  clusters  by  clicking  on  the  computer 
screen  at  the  centers  of  the  clusters.  These  2-D  models  will  be 
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refined  by  EM  algorithm  and  compete  according  to  MDL 
criterion  or  human  justification.  The  winning  model  in  2-D 
space  will  be  transferred  back  to  original  data  space  to 
initialize  the  data  model  in  that  space.  Then  EM  algorithm  in 
original  data  space  will  refine  the  model  and  obtain  the 
partition  of  data  at  that  level.  When  no  more  new  clusters  can 
be  found  in  the  model  validation  step,  the  algorithm  ends  and 
a  final  partition  is  obtained. 

III.  Results 

A.  Data  treatment 

We  applied  nICA  to  Saccharomyces  cerevisiae  gene 
expression  dataset  [9].  The  dataset  contains  6152  genes  with 
Open  Reading  Frames  (ORFs)  and  173  samples  that  include 
the  different  experimental  conditions  [9] :  temperature  shocks, 
amino  acid  starvation,  and  progression  into  stationary  phase 
etc.  As  in  [3],  we  also  used  KNNimpute  to  fill  in  missing 
values.  And  due  to  the  triviality  of  clustering  environmental 
stress  response  (ESR)  genes  defined  by  [9],  we  eliminated 
them  in  our  analysis.  The  final  dataset  contains  5284  genes 
and  173  samples.  To  evaluate  the  experimental  results,  we 
measured  the  biological  significance  of  each  cluster  using 
Gene  Ontology  (GO)  annotation  database.  We  mainly 
measured  each  cluster  of  nICA  in  terms  of  the  biological 
process  and  molecular  functional  categories  using  /7-value 

[3]. 
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Fig.  2  Comparison  of  the  scatter  plots  of  the  first  three  independent 
components  from  nICA/ICA.  (a)  Results  from  nICA;  (b)  Results  from 
ICA.  Each  sub-panel  shows  the  two  subsequent  components  plotted 
against  each  other.  In  (a),  process-specific  genes  are  highly  biased  on  two 
non-negative  axes,  whereas  the  results  of  ICA  in  (b)  are  not. 


B.  Experiment  and  results  evaluation 

To  use  nICA  model,  we  took  the  inverse-logarithm  of  the 
data  before  further  analysis.  Hence,  in  our  nICA  model,  the 
microarray  expression  corresponds  to  a  linear  additive  model 
of  interactions  among  biological  processes.  Since  our  goal  is 
to  find  some  most  relevant  components,  dimension  reduction 
using  PCA  was  first  applied  to  the  data  with  90%  of  energy 
maintained.  Then  nICA  was  applied  to  the  dimension  reduced 
dataset.  As  a  result,  we  obtained  the  eighteen  non-negative 
independent  components.  For  comparison,  we  also  did  the 
experiment  using  an  ICA  algorithm  with  the  same 
parameters. 

Fig.  2  shows  the  first  3  independent  components  from 
nICA  and  ICA  respectively.  It  is  clear  that,  comparing  to 
ICA,  nICA  is  effective  in  separating  sources  as  independent 
non-negative  “biological  processes”  in  which  process-specific 
genes  are  highly  biased  onto  two  orthogonal  axes  respectively 
showed  in  each  sub-panel  in  Fig.  2. 

In  Table  1  we  only  listed  seven  most  significant  clusters 
resulted  from  our  nICA  and  VISDA  approach.  We  measured 
the  biological  significance  of  each  cluster  using  GO 
annotation  database.  The  p-\ alue  of  each  cluster  was 
calculated  according  to  its  overlap  with  the  functional 
annotations  in  GO  (see  [3]  for  the  details).  Among  those 
functional  categories  detected  significantly  by  both  nICA  and 
ICA  clusters,  there  are  five  out  seven  clusters  that  nICA 
produced  significant  lower  /7-values  than  ICA  did.  From  these 
experiments,  it  seems  to  us  that  nICA  followed  by  VISDA 


can  extract  more  coherent  groups  of  genes  in  terms  of  their 
functional  categories. 

To  further  evaluate  our  nICA-based  clustering  method,  we 
used  the  z-score  introduced  in  [8]  to  conduct  a  comparative 
study.  As  described  in  [8],  the  z-core  is  based  on  the  mutual 
information  between  clustering  results  and  the  gene 
annotation.  The  higher  scores  indicate  clustering  results  more 
significantly.  We  compared  the  clustering  results  of  nICA  and 
ICA  under  the  same  parameters  and  the  z  scores  are  shown  in 
Fig.  3.  As  we  can  see  from  Fig.  3,  nICA  algorithm 
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Fig.  3  Z  score  for  the  nICA  (asterisk)  and  ICA  (circle).  At  each  level  of 
the  hierarchy  of  VISDA,  we  recorded  all  the  intermediate  clusters.  At 
last,  we  got  41  clusters  for  nICA  and  28  clusters  for  ICA.  We  inputted  all 
these  clusters  to  compute  the  z  scores  and  drew  the  curves  here.  So  based 
on  the  figures,  it  is  reasonable  to  draw  the  conclusion  that  clustering 
methods  by  nICA  has  found  a  finer  structure  than  ICA  has. 
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TABLE  I 

The  seven  most  significant  clusters  of  NON-NEGATIVE  ICA 


Cluster  index  Gene  Ontology  term 

Cluster  frequency 

Genome  frequency  of  use  P-value 

6 

Ty  element  transposition 

(24  /  73,  32.8%) 

(95/7291,  1.3%) 

3.7E-27 

DNA  transposition 

(24  /  73,  32.8%) 

(108/7291,  1.4%) 

7.3E-26 

DNA  recombination 

(24/73,  32.8%) 

(192/7291,2.6%) 

4.2E-20 

RNA-directed  DNA  polymerase  activity 

(14/73,  19.1%) 

(52/7291,0.7%) 

2.2E-16 

DNA-directed  DNA  polymerase  activity 

(15/73,20.5%) 

(67/7291,0.9%) 

2.5E-16 

16 

glycolysis 

(9/109,8.2%) 

(21  /7291,0.2%) 

4.5E-11 

Glucose  metabolism 

(12/109,  11.0%) 

(65/7291,0.8%) 

3.6E-10 

17 

proteolysis 

(40/250,  16%) 

(164/7291,2.2%) 

4.3E-22 

ubiquitin-dependent  protein  catabolism 

(24/250,  13.6%) 

(128/7291,  1.7%) 

5.5E-20 

endopeptidase  activity 

(25/250,  10%) 

(62/7291,0.8%) 

4.5E-19 

22 

amino  acid  and  derivative/metabolism 

(21  /43,48.8%) 

(199/7291,2.7%) 

8.5E-22 

amino  acid  metabolism 

(20/43,46.5%) 

(183  /7291,2.5%) 

5.4E-21 

oxidoreductase 

(10/43,23.2%) 

(247/7291,3.3%) 

1.4E-06 

23 

amino  acid  biosynthesis 

(19/85,22.3%) 

(102/7291,  1.3%) 

1.0E-17 

catalytic  activity 

(43  /  85,  50.5%) 

(1937/7291,26.5%) 

2.1E-06 

24 

cellular  response  to  nitrogen  starvation 

(4/47,8.5%) 

(5/7291,0.0%) 

3.9E-08 

cellular  response  to  nitrogen  levels 

(4/47,8.5%) 

(5/7291,0.0%) 

3.9E-08 

asparagine 

(4  /  47,  8.5%) 

(5/7291,0.0%) 

3.9E-08 

33 

generation  of  precursor  metabolites  and  energy 

(32  /  68,  47.0%) 

(231  /7291,3.1%) 

8.8E-30 

oxidative  phosphorylation 

(19/68,27.9%) 

(46/7291,0.6%) 

4.0E-26 

hydrogen  ion  transporter  activity 

(20  /  68,  29.4%) 

(55/7291,0.7%) 

2.1E-26 

The  selected  clusters  are  listed  along  with  the  functional  categories  with  the  smallest  p-value.  Numbers  in  parentheses  in  the  third  column  show  the  number 
and  percentage  of  genes  within  the  cluster  that  are  presented  in  one  of  the  functional  category.  For  instance,  (24/73,  32.8%)  means  the  cluster  has  73  genes, 
among  which  24  (32.8%)  genes  are  annotated  with  “Ty  element  transposition”.  And  the  numbers  in  the  fourth  column  are  presented  in  the  similar  way 
which  corresponds  to  the  total  number  within  the  whole  genome  set  that  are  annotated  with  one  of  the  special  categories  in  GO  system. 

consistently  performed  better  than  ICA  with  an  average 

increase  of  z-score  of  5.  References 


IV.  CONCLUSION  AND  DISCUSSION 

This  paper  presents  a  new  gene  clustering  approach, 
namely  nICA-based  approach  for  composite  module 
discovery.  By  projecting  the  gene  expression  data  onto  nICA 
space,  co-regulation  structure  of  modules  are  revealed  and 
highlighted.  Using  a  pre-screening  and  VISDA  clustering 
procedure,  we  can  identify  biological  process  enriched 
clusters  with  coherent  functional  annotations.  The 
experimental  results  on  a  yeast  data  set  have  demonstrated  its 
advantages  over  conventional  ICA-based  approach. 

Although  nICA-based  approach  exhibits  some  promise 
for  gene  clustering,  there  is  future  work  to  be  conducted.  For 
example,  we  notice  that  in  the  de-mixing  matrix  W,  there  are 
some  negative  values  that  need  to  be  properly  explained,  i.e., 
how  these  composite  modules  are  involved  in  the 
corresponding  biological  process.  Another  possible  direction 
is  that  we  may  perform  gene  clustering  to  find  groups  of 
genes  under  distinctive  regulators  or  combinations  of  genes 
of  these  regulators. 
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